本発明者らは、abinitio分子動力学を用い、かつ新たな変分計算法を開発し、(i)各プロセス要素の揺らぎや、各要素のずれの重畳現象も診断し、(ii)所望のシーケンス通り実行されているかのチェックや、(iii)変動因子の発見も可能とするシミュレーションシステムを開発した。かかるシステムのさらなる特徴は、(vi)プロセスの変動を数学的・統計的に、しかも効率よく取り扱えるようにした点と、(v)新しい推論エンジンを付加した点と、(vi)新しい逆方向システムを付加した点にある。
即ち、本発明は、複数の工程からなる半導体装置の製造方法において、前記複数の工程の少なくとも1つにおける製造装置内に形成された膜、半導体基板上に形成された膜、半導体基板表面、及び半導体基板内部の少なくとも1つに関する複数の特性の少なくとも1つを測定して実観測データを得る工程と、abinitio分子動力学プロセスシミュレータ又は経験的ポテンシャルを与えた分子動力学シミュレータにより、前記複数の工程の少なくとも1つにおける、複数の製造工程因子の少なくとも1つの設定値、及びこの設定値のゆらぎを入力値とし、この入力値を基に変分計算を行い、前記複数の特性の少なくとも1つの予測データを得る工程と、前記予測データと実観測データとを逐次、実時間で比較検定する工程と、前記比較検定により、前記複数の製造工程因子の少なくとも1つの設定値と、前記実観測データから推測される前記複数の製造工程因子の少なくとも1つとの間に有意差が認められた場合、前記製造工程因子を逐次実時間で修正処理する工程とを具備することを特徴とする半導体装置の製造方法を提供する。
本発明は、上述の半導体装置の製造方法(請求項1)において、前記abinitio分子動力学プロセスシミュレータ又は経験的ポテンシャルを与えた分子動力学シミュレータは、前記複数の工程において形成又は処理される半導体材料、絶縁材料、又は導電性材料を構成する原子の運動が平衡に達したところで、時刻を0に設定し、その後、時刻tまでをサンプリング時間とし、各原子の分極率、電荷、及び位置ベクトルとの時間相関を求め、所定時間における光学的スペクトルを求める機能を有し、前記abinitio分子動力学プロセスシミュレータ又は経験的ポテンシャルを与えた分子動力学シミュレータから導かれる予測データは、個々の原子の時々刻々の運動から導かれる振動挙動、その固有振動数、及び固有振動数の分布情報を含むことを特徴とする半導体装置の製造方法を提供する。
本発明は、上述の半導体装置の製造方法(請求項1)において、前記abinitio分子動力学プロセスシミュレータ又は経験的ポテンシャルを与えた分子動力学シミュレータは、前記複数の工程において形成又は処理される半導体材料、絶縁材料、又は導電性材料を構成する原子の運動から、該材料の材料強度定数を逐次算出し、これを用いて該材料に誘起される応力を求める機能を有し、該応力が規定値を越える場合には、前記製造工程因子の少なくとも1つが逐次修正処理されることを特徴とする半導体装置の製造方法を提供する。
本発明は、上述の半導体装置の製造方法において、前記半導体材料、絶縁材料、又は導電性材料を構成する原子は、シリコン原子及び酸素原子である半導体装置の製造方法を提供する。
本発明は、複数の工程を実施する半導体装置の製造装置において、前記複数の工程の少なくとも1つにおける製造装置内に形成された膜、半導体基板上に形成された膜、半導体基板表面、及び半導体基板内部の少なくとも1つに関する複数の特性の少なくとも1つを測定して実観測データを得る手段と、abinitio分子動力学プロセスシミュレータ又は経験的ポテンシャルを与えた分子動力学シミュレータにより、前記複数の工程の少なくとも1つにおける、複数の製造工程因子の少なくとも1つの設定値、及びこの設定値のゆらぎを入力値とし、この入力値を基に変分計算を行い、前記複数の特性の少なくとも1つの予測データを得る手段と、前記予測データと実観測データとを逐次、実時間で比較検定する手段と、前記比較検定により、前記複数の製造工程因子の少なくとも1つの設定値と、前記実観測データから推測される前記複数の製造工程因子の少なくとも1つとの間に有意差が認められた場合、前記製造工程因子を逐次実時間で修正処理する手段とを具備することを特徴とする半導体装置の製造装置を提供する。
本発明は、abinitio分子動力学プロセスシミュレータ又は経験的ポテンシャルを与えた分子動力学シミュレータにより、半導体装置の製造プロセスの複数の工程の少なくとも1つにおける、複数の製造工程因子の少なくとも1つの設定値、及びこの設定値のゆらぎを入力値とし、この入力値を基に変分計算を行い、前記複数の工程の少なくとも1つにおける製造装置内に形成された膜、半導体基板上に形成された膜、半導体基板表面、及び半導体基板内部の少なくとも1つに関する複数の特性の少なくとも1つの予測データを得る工程と、この予測データと前記複数の工程の少なくとも1つにおける前記複数の特性の少なくとも1つを測定して得られた実観測データとを逐次、実時間で比較検定して、前記製造工程因子の診断及び判定を行う工程とを具備するシミュレーション方法を提供する。
本発明は、上述のシミュレーション方法(請求項6)において、前記比較検定により前記複数の製造工程因子の少なくとも1つの設定値と、前記実観測データから推測される前記複数の製造工程因子の少なくとも1つの間に有意差が認められた場合、前記製造工程因子の修正値を実時間で求める工程を具備することを特徴とするシミュレーション方法を提供する。
本発明は、半導体装置の製造プロセスの複数の工程の少なくとも1つにおける製造装置内、半導体基板表面、及び/又は半導体基板上に形成された膜に関する複数の特性の少なくとも1つを測定して実観測データを得る手段と、abinitio分子動力学プロセスシミュレータ又は経験的ポテンシャルを与えた分子動力学ミュレータにより、前記複数の工程の少なくとも1つにおける、複数の製造工程因子の少なくとも1つの設定値、及びこの設定値のゆらぎを入力値とし、この入力値を基に変分計算を行い、前記複数の特性の少なくとも1つの予測データを得る手段と、前記予測データと実観測データとを逐次、実時間で比較検定する手段とを具備する、半導体装置の製造工程因子を診断・判定するシミュレータを提供する。
本発明は、abinitio分子動力学プロセスシミュレータ又は経験的ポテンシャルを与えた分子動力学シミュレータにより、半導体装置の製造プロセスの複数の工程の少なくとも1つにおける、複数の製造工程因子の少なくとも1つの設定値、及びこの設定値のゆらぎを入力値とし、この入力値を基に変分計算を行い、前記複数の工程の少なくとも1つにおける製造装置内に形成された膜、半導体基板上に形成された膜、半導体基板表面、及び半導体基板内部の少なくとも1つに関する複数の特性の少なくとも1つの予測データを得る手段と、この予測データと前記複数の工程の少なくとも1つにおける前記複数の特性の少なくとも1つを測定して得られた実観測データとを逐次、実時間で比較検定して、前記製造工程因子の診断及び判定を行う手段とを具備する、半導体装置の製造工程因子を診断・判定するシミュレータを提供する。
一方、本発明者等は分子動力学シミュレータについて、前記したような非常に計算が遅いという問題、及び対象とする物質の実行領域が、非常に小さいという問題に対してそれぞれ新しい考えを導入して克服し、分子動力学シミュレータとプロセス変動算出プログラムとを新しい方法で合体させ、しかも、この合体により、今までにない新しい機能を引き出すことが出来た。以下に順に説明する。
まず(1)の計算速度の問題に関して少し詳しく触れてみる。分子動力学では、原子一個一個について、いわばニュートン力学に従って、運動方程式を立てて、これを、着目している100個または1000個の原子全部に渡って連立させて解くものである。
この運動方程式の基になる力の計算が一番複雑である。これが、従来の計算手法では遅かったゆえんである。各原子は、その廻りの原子から力を受けている。
近くの原子からも、遠くの原子からも、力を貰っている。これを考慮する必要がある。一般にはテレソフの式に従って導出されるポテンシャルエネルギVijをもとめ、これを基に、力を算出して行く。しかし、従来では、概ね、ポテンシャルを、距離や方向等をもとに、予め各原子に加わる力について方向毎のテーブルを作っておく手法を取っていた。この方法では、上記の100個または1000個の原子について、毎回、テーブルを参照し、ポテンシャルの大きさを見て行くものである。これは、非常に時間の掛かるものである。またポテンシャルそのものからは、力を算出することは出来なくて、ポテンシャルの大きさのそれぞれの方向に対する傾きが、それぞれの方向の力に成るわけである(図2参照)。これゆえ、非常に計算機の時間がかかることが分かる。
なぜこの様に、数値テーブルを作るかと言うと、ポテンシャルの形状を表した式が非常に複雑であるからである。上にも述べた様に、ポテンシャル(V)は、廻りの原子の位置はお互いの角度(θ)、さらには、距離(r)等が複雑にからみ合っている。これ故に、これを、数学的に、微分して関数を作るのは非常に難しい。例えばθはrの関数であるのに、また廻り回って再びrの関数に成ったりする。
本発明者らは、この様な一般に非常に複雑で、上記の様にポテンシャルの表現式の引数が、巡回する様な場合について、従来、偏微分が試みられなかった部分において、新しい高階巡回偏微分式を作成した。これにより、式自体の変形は一見非常に複雑になるが、計算は非常にはやく出来るようになった。従来のように、数表を参照しなくて良いようになった。
では、この高階巡回偏微分式の一番の概念を以下に示す。これは、一般に3個の原子を想定して、考えている。すなわち、原子3個に関して一般的に記述すれば、あとは、この手法を広げて行けばよい。3個の原子があれば、各原子との距離や、その着目している原子の間の角度等も算出できる。これらにより、ポテンシャルを一義的に求められ、これをもとに、何れの方向にも、ポテンシャルの「傾き」を容易に算出することが出来る。本発明者等は、巡回して用いている変数について全て偏微分式を作っておき、これらを、連立させて解くことで、正しく、微分型が作られていることを確認した。
以下本発明によるシステムの概要を具体的に説明する。最初に、図1を参照して、本発明の応用の実際を説明する。即ち、プロセス装置Vは、実際に稼働している半導体製造システムの一部である。ここでは、スペクトルメータVIによって光の吸収スペクトルをみている。
スペクトルメータVIによって得られた光の吸収スペクトルは、一旦ファイル4として格納され、スケジューラVIIに送られる。スケジューラVIIには、後述する本発明による予測値も入力され、これらは対比診断検定部IVで比較される。
対比診断検定部IVでの比較結果が、許容値を超えていれば、スペックを満たさないものとしてラインから外される。また、満たしていれば、その変動分を推論エンジンとしてのプロセス変動算出ユニットIに送られる。
プロセス変動算出ユニット1では、基本的には原子レベルまで考慮しないマクロ的な性質を計算し、ファイル2に格納する。ファイル2の内容は、分子動力学シミュレータIIで原子間位置ベクトルと結合角が算出される。
分子動力学シミュレータIIで算出された原子間位置ベクトルと結合角は、膜特性算出ユニットIIIに入力され、ここで基本膜特性が算出されてファイル3に戻される。膜特性算出ユニットIIIで算出された基本膜特性は、プロセス変動算出ユニットIにも戻され、パラメータの一部とされる。
このように、本発明によって、半導体装置製造プロセスを、テストピースなしに、所望の工程どおり、または修正しながら進行する事を可能とする半導体装置の製造方法が実現される。
実際には、分子動力学シミュレータIIの内容が特に重要である。具体的には、図7のブロック図に示したように、2つの部分からなっている。第1の部分IIaは、本件発明者によって導入された高階巡回偏微分法を利用するものであり、第2の部分IIbは、三次以上の高次項を考慮した解析を行うものである。なお、以後、高階巡回偏微分法や三次以上の高次項を考慮した解析の詳細が順次説明される。
即ち、本発明者らは、純粋数学理論にたち帰り、新たな変分計算法を確立するとともに、新しい計算化学理論により、各種スペクトルをも予測する技術を開発した。以下、本発明における新しいabinitio分子運動学について説明する。
本発明者等が開発した新しいシミュレーションの全体構成を図3に示す。分子動力学部分(MD)と密度汎関数(DFT)とを合成した新しい手法を提案した。まず本発明者らによる分子動力学法と密度汎関数法との使い分けを図3を用いながら以下に示す。
図3に示すように、分子動力学(MD)部分はフローチャートの下半分を占め、密度汎関数部分(DFT)は上半分を占めている。分子動力学法は零Kでない所謂有限温度でのやや大きな系を扱い、原子(イオン)間のポテンシャルには、予め第一原理から求めたポテンシャルを用いる。他方、密度汎関数(DFT)は、小さな系の基底状態を求める時に使うことにした。そして、当然DFTの部分は温度は零Kとした。ポテンシャル部分は、電子系の基底状態から求めた。具体的には波動関数(KS方程式)を解き、局所密度汎関数を併用した。
DFT部分で、縮重のない電子系の基底状態の全エネルギーを求めて置く。全エネルギーは(波動関数までさかのぼらずに)電子密度n(r)の汎関数E[n(r)]である。
即ち、E[n(r)]は、下記式(1201)により表される。
右辺第1項のVz (r)は核のクーロン場などの外場、第2項は電子間のクーロン相互作用エネルギー、第3項は、電子の1体近似を想定して、電子間相互作用のない時の運動エネルギーTs[n(r)]であり、下記式(1202)により表される。
第4項は交換相関エネルギーExc[n(r)]で、その形は局所密度近似による。電子密度n(r)は下記式(1203)により表される。
基底状態は、式(1201)のE[n(r)]の変分を0とおいて得る、下記式(21),(22)に示すKS方程式(Kohn−sham equation)を解いて、固有値Eの小さい方からN個採用する。式(1211)にn(r)が入っているから、n(r)、ψ(r)、有効ポテンシャルVeff (r)を自己無撞着(self−consistant)に解かなければならない。
〔数4〕
[−∇2+Veff(r)]Ψ(r)=EΨ(r) …(1211)
Veff(r)=VZ(r)+∫2n(r′)/(|r−r′|)・
dr′+δExc[n(r)]/δn(r)
…(1212)
交換相関エネルギーExc[n(r)]は、局所密度近似(LDA)で1電子当たりのエネルギーεxc(n)に電子密度nを掛けて積分することにより、下記式(1213)を得た。
〔数5〕
Exc[n(r) ]=∫εxc(n(r))n(r)dr …(1213)
εxc(n)の形は、Wigner他の式など、いろいろ提案されている。更に、内殻電子は考えずに価電子だけ扱う時には、式(1211)のVz (r)として核のクーロン場をそのまま使うのでなく、核の位置での特異性を消したノルム保存擬ポテンシャルでおきかえる。このようにならすことによって、フーリエ展開での成分の個数をおさえることができた。
全エネルギーとして、本発明者等は厳密を期するため、原子(イオン)間のクーロン・ポテンシャルも含むものを使った。即ち、下記式(1221)である。
ここで{RI}はイオンの核の座標、{αν}は外部の制約条件(体積Ω、ひずみεLνなど)。式(1221)の右辺第1項は、DFTの式(1201)から右辺のTs [n(r)]であり、式(1201)右辺の残りが式(1221)のUの1部である。その部分は、DFTの場合と同様に1体近似の式(1211)からVeff(r)で表され、やはりVz (r)を擬ポテンシャルで置き換えて使う。
ここで新しい考え方を入れた。即ち、Eをポテンシャル・エネルギーと見なして、他に、下記式(1222)に示す仮想的な運動エネルギーを導入する。
従ってラグランジアンは、下記式(1223)により表される。
L=K−E …(1223)
式(1222)中のM
Iは、本当の原子(イオン)の質量だがμ、μνは仮想的な質量であり、後述のように目的に応じて変えた。束縛条件として、下記式(1224)に示す正規直交条件
〔数8〕
∫ψ
i *(r,t)ψ
k(r,t)dr=δ
ik …(1224)
を課するから、オイラー方程式を作る時に未定乗数Λ
ikを導入して、Lに下記式(1231)を加えたものの変分をとる。その結果、下記式(1232),(1233),(1234)が得られる。
式(1222)の運動エネルギーの温度Tが対応する。古典統計力学のエネルギー等分配則によって各自由度が温度を持っているのだから、Σ1/2 MI RI 2 との関係でいえば、仮想的でなく現実の温度である。またμ、μνの大小によってψi 、ανを抑制したり増強したりすることが出来る。
例えば、μ<<MIとして高温からT→0とすれば、RI をとめたままψi が変化して電子系の基底状態を得る。この時、式(1231)の左辺が0になり、右辺は式(1221)とその下の注意により、下記式(30)となるから、ψi の線形結合を適当に作る{Λik}が対角線化されてKS方程式(1235)を得る。
〔数10〕
[∇2−Veff(r)]Ψj(r)+ΣΨk(r) …(1235)
つまりDFTでKS方程式を解いたのと同じ結果をsimulated annealingで得たことになる。ただし、温度Tの状態を作るのに、モンテカルロ法によらず、仮想的な力学系の運動を追っているから、dynamical simulated annealing(DSA)と呼ばれる。
式(1231)の積分をする時、必ずしもψi をそのまま変数として扱うのでなく、それを平面波で展開した展開係数(つまりフーリエ係数)を変数として扱うことができる。この展開係数の個数をむやみにふやさないためには、DFTの場合と同様に、擬ポテンシャルで核の位置での特異性を消しておくことが必要である。特に、結晶の周期性を利用できるときには、上記の{Λik}を対角線化するψi の線形結合をφnkと記す。kはブリユアン帯域内の波数ベクトル、nは帯域インデクスである。
φnkは、下記式(1241)を満たすはずである。
〔数11〕
μS
nk(r,t)=[∇
2−V
eff(r)]φ
nk(r)+λ
nkφ
nk(r,t)
…(1241)
λnkは行列{Λ
ik}の固有値であり、エネルギー準位である。φnk(r,t)を展開して、下記式(1242)が得られる。
擬ポテンシャルが局所的であるならば(この条件は一般に満たされない)、有効ポテンシャルV
effも展開して、下記式(1243)が得られる。
式(1242)と式(1243)を式(1241)に代入して、下記式(1251)が得られる。
右辺の最後の頃でK=K′+K″とおくと、Σexp [i(K+G)r]が全部の項に共通となるから外して、下記式(1252)を得る。
〔数15〕
μTn,K+G(t)=[−|K+G|2+λnk]Cn,K+G(t)
−ΣVG-G′(t) Cn,K+G′(t)
=−[|K+G|2+V0(t)−λnx]Cn,K+G(t)
−ΣVG-G′(t) Cn,K+G′(t)
…(1252)
ここで周波数を、下記式(1253)に示すように定義し、
〔数16〕
ω=((|K+G|2+V0(t)−λnx)/μ)1/2
…(1253)
式(1251)の左辺を差分になおすと、下記式(1261)が得られる。
〔数17〕
Cn,K+G(t+Δt)=2cos(ωΔt)Cn,K+G(t)−Cn,K+G(t−Δt)
−2[1−cos(ωΔt)]{ΣG´VG-G´(t)・Cn,K+G´(t)}
…(1261)
この式は、振動部分の積分を解析的にすましているので、純粋に数値的な方法よりもΔtを大きくすることにより、計算量を減少させることが出来る。
クーロンポテンシャルは、以下のようになる。即ち、本発明者は、クーロンポテンシャルを分解してゆくと4項に分かれる事を見い出した。特に、従来は3項しかなかったが、新しく第4項があることを確認した。これらの方法を順次説明する。
本発明者らは、まず厳密な式の出発として、誘電率を加味した基本式から解き始めた。まず基本式である下記式(1262)から解き始めた。
導体ε=∞と真空ε=1の場合とでクーロンポテンシャルが異なり、Lは(立方体の)単位結晶の一稜、Σは単位結晶内でとり、Zi 、ri は、それぞれ第i粒子の荷電と位置である。これは球内の荷電によって球の内面に分極が生じることによる。導体でない球の内面に双極子の層が出来るが、上記式の最終項がちょうどその効果を打ち消す働きをする。Ewaldの方法は左辺、ε=∞のものを与えるから、真空内の値を求めるには上記式の最終項を加えなければいけない。
結果だけ掲げる。
クーロン・ポテンシャルを下記式(1271)により表わすものとする
式(1271)において、Nは単位結晶内の原子数、その中の第1原子の位置と荷電がri ,Zi であり、nは単位結晶とそれを周期的にずらしたものを指定するベクトルであって、下記式(1272)に示すように設定した。
n=nxZx+nyZy+nzZz …(1272)
と表される。ここに、Zx ,Zy ,Zz がそれぞれx,y,z方向の稜のベクトルで、nx ,ny ,nz はそれぞれ(バルク結晶の場合)−∞から+∞にわたる整数である。Σ′の′はn=0の時のj=iを除くことを意味するものとする。
ここで本発明者らは下記式(1273)に示す新しいF関数を導入する。
上記式(1273)において、G(r,t)はtの周期関数である。これはフーリエ級数で表現できることを見い出した。
G(r,t)は、更に下記(1281)に示すように変形することが出来る。
上記式(1281)において、mは逆格子ベクトルであり、下記式(1282)により表わされる。
〔数22〕
m=m
x・(1/L
x,0,0)+m
y・(0,1/L
y,0)
+m
z・(0,0,1/L
z)
…(1282)
指数mx ,my ,mz はそれぞれ−∞から∞にわたる整数である。m=0の時は、exp […]=1だから、ΣZi =0であり、単位結晶内の総荷電が0であれば、m=0の項は消える。あらかじめm=0の項を除くと、下記式(1283)となる。
本発明者等は積分範囲をkで分け、G(r,t)の2つの形を使い分けて、下記式(1291)を得る。
ここで本発明者等は、公式として下記式(1292)を用いた。
これにより、始めのクーロンポテンシャルは、下記式(1293)により表わされる。
Σ′の右上の′印は、n=0の時、j=iを除くことを意味する。
この式(1301)はrを含まないから、力には関係ない。この式(1301)は、更に下記式(1302)に示すように変形することが出来る。
以上の結果をまとめると、クーロンポテンシャルは、結局、下記式(1311)〜(1315)により表わされる。
ここに、nは、n=nx ・(Lx ,0,0)+ny ・(0,Ly ,0)+nz ・(0,0,Lz )であり、、mは、m=mx ・(1/Lx ,0,0)+my ・(0,1/Ly ,0)+mz ・(0,0,1/Lz )である。上記の表現の具合のよい点は、もとのΣの項が逆数のオーダーでしか減衰しないのに対して、Φ(1) ではΣの項がerfcの因子によって、Φ(2) ではΣの項がexp の因子によって、急速に減衰することである。Φ(3) での減衰との遅速に逆向きに効くから、両者のバランスがとれるような適当なκを選ぶ必要がある。これらは、距離の近いところから順にクーロン力の寄与を計算した結果であり、しかも周囲が導体である場合の結果である。
周囲が真空である場合によればもう1項が加わる。これが特に従来省みられなかった部分である。
一般に分子動力学では、原子間ポテンシャルが最も大切な量である。これには、現状では、単純には、2体間でのみ記述するものも結構多い。これは、分子動力学計算が非常に、時間のかかるもので、その時間節減のために、ポテンシャルを簡便なものを使っているわけである。この様に、ポテンシャルを簡便にしてしまうと、計算自体は速くなるが、実際を反映したものではありえない。
上に記したように重要な原子間ポテンシャルは、一般には、3体で記述できる。これの意味合いは以下の通りである。即ち、第一の原子(i)と、第二の原子(j)との作用に、さらに第三の原子(k)から第二(j)を通して入ってくる効果がある。これを図示したものが図4である。図に有るように、rを粒子間の距離とし、θを粒子角度とする。これらの諸量の主従の関係は、図5に示す通りである。特に図5で分かるように、非常に複雑であることが分かる。一般に、ポテンシャルを微分したものが、力になるわけであるが、これで分かるように計算は、非常に複雑になる。ここでもう一度、従来の計算手法と本発明の位置つけを図6を用いて示した。
本発明者らは、図6にある高階巡回偏微分法を新しく開発して、これにより、高精度に且つ高速に演算できる様にした。本発明者等の方法がどれくらい速いかを以下に実際の比較例を入れて記す。まず、図6にある簡便な2体間ポテンシャル計算に関しては、ここでは、触れない。なぜなら、半導体LSIの物理現象を記述するには殆ど使用に耐える精度が得られないので、ここでは、議論しないことにした。一方、3体問題を入れたポテンシャルで、従来の方法として、テーブル参照方法のプログラムを、本発明者等は作成した。これによると、rijと,θijkと、rjkに関して、表を作成した。実際には、block dataとして、配列にして格納した。そして、この参照項目からずれていた時のことを考え、補完プログラムをも作成した。補完には、cubic spline法を用いた。ある値、即ち、rijと, θijkと、rjkとに関して、テーブルを引用出来るようなプログラムにした。し かし、所望の条件(rijと,θijkと、rjk)が、配列内のどの行列に入るかを探る必要がある。本発明者等は、出来るだけ高速に成るように気を配りながらプログラムを作成したが、やはり、数字の大きさによる場合分けをする部分が必要で、そこにはif then else文を用いた。この様にして、所望の条件(rijと,θijkと、rjk)がテーブルのどの枠にはまるかを、見いだすとともに、その中での補完をおも行った。これで、まず、所望の条件(rijと,θijkと、rjk)でのポテンシャルが求まったわけであるが、実は欲しいのは、この点における力である。
本発明者等は、微小変化させた部分でもう一点についても、ポテンシャルを同様の手続きを用いて求め、両者の差から、平均変化率の考えを用いて、微分、即ち力を求めた。
これで分かるように、非常に多くのif then else手続きを用いた。電子計算機では、if then else手続は比較的苦手で、遅い。事実、本発明者等が作成した、参照型プログラムを用いた場合、力を計算する部分のみについて、詳細について調べた。その結果、1600個のSi原子からなる単結晶基板で、絶対温度300Kで、外圧1気圧の場合であるが、参照型プログラムでは、約350倍の時間を要した。これで分かるように、本発明者等が提案する、高階巡回偏微分法のほうが、有利であることがわかる。また参照型プログラムでは、原子が異なる度に、その都度データ表を作成する必要がある。
次に分子動力学から材料特性値等をどの様に抽出して、それを汎用シミュレータにどの様に入れているのか」に関しても簡単に説明する。
通常、LSIプロセスに於いては、種々様々な評価技術を使っている。赤外吸収による膜の物性値評価法、さらには、ラマン分光による膜中の応力評価、X線による膜の構造解析等々がある。またこれら光学的なものに限らず、膜の誘電率等電気的な評価技術もある。
本発明者等は、それぞれの評価手法の原理を、一旦、鋭意、理論的に個々原子の運動にまで分解し詳細に追跡した。その上で、本発明者等の開発した分子動力学を用いて、個々原子の動きを追跡し、それぞれの評価量を初めて理論的に算出し、たとえばスペクトルの様な形に直接表示することが出来た。これにより、例えば、実測値との定量的な比較が出来る等、従来にない大きな技術の進捗を得たわけである。
本発明者等が、初めて克服した理論的手段による原子運動の追跡と、具体的な測定評価量の導出までを以下に記す。幾つかのものがあるが、基本的には、本発明者が導出した基本概念は以下の通りである。即ち、着目する原子のモーメントを計算する。そして、これを、着目する時間の窓の中で積算し、しかる後に、フーリエ変換し、固有スペクトルを得るものである。これが大きな骨子である。
次に、Siを例に、経験的なポテンシャルに関して、本発明者がどの様にして、そのポテンシャルを力に変換したかを示す。これは、一例であるが、Siに限らず、イオン結合による効果を考慮しなければ、例えば酸化膜のポテンシャルも、このSiの問題と同様に、3体問題として、ここに述べる高階巡回偏微分を用いることが出来る。本発明者等は、Siに関してどんなポテンシャルが良いかも吟味した。
ここで用いるポテンシャルはTersoffによるものである。(Phys.Rev.Lett.56.632 (1986), Phys.Rev.B37,6991(1988))。
今、Tersoffによれば、i番目のSiに関する全ポテンシャルは、
で記述できる。これは3体で記述しているので、本発明者等が注意を喚起したいのは、上記式(0001)に於いて、Vij≠Vjiである。着目するSi原子の位置番号をiとし、その周辺の他の粒子番号をjとすると、上記Vijは、
〔数31〕
Vij=fc(rij)[aijfR(rij)+bijfA(rij)] …(0002)
である。ここにrは粒子間の距離である。また、fc(rij)は、カットオフ関 数で、fR(rij)は斥力を示し、fA(rij)は引力を示す。aijは配位数を考慮したカットオフ係数、bijも配位数を考慮したカットオフ係数である。配位数を表現するパラメータを用いて3体的表現を行っていることになる。これは、他のものでも基本形は同じであるので、ここでは、本発明者等が初めて開発した巡回手法の思想を交えながら以下に示す。
fR(rij)とfA(rij)は、Morse型のポテンシャルを変形したもので、fR(rij)=Aexp(-λ1r)、fA(rij)=−Bexp(-λ2r)である。
この内、λ1とλ2は、定数であり、その大きさは、原子間距離程度の値の逆数である。
ところで、カットオフ関数fc(rij)とは、
〔数32〕
Vij=fc(rij)[aijAexp(−λ1rij)−bijBexp(−λ2rij)]
fc(r)=1 (r≦R−D)
fc(r)=1/2−1/2sin[π/2(r−R)/D](R−D<r<R+D)
fc(r)=0 (r≧R+D)
…(0003)
であり、ここに、Rは通常対象とする構造の第一隣接ゾーンだけを含む様にその寸法を選ぶ。その値は、本発明者の詳細な吟味によると大体2〜3Åである。次に、bijは実配位数を示すことになるが、ここでも上記カットオフ関数を使う。
その定義は、Tersoffによれば、 bij=(1+βnζij n)-1/2n …(0004)
である。ここに、
〔数33〕
ζij=Σfc(rik)g(θijk)exp[λ3 3(rij−rik)3] …(0005)
である。Σ記号はk≠i,jで回す。ここで分かるように、ζijの意味は第3の原子kが入ることによる環境因子であるので、iから見た場合と、jから見た場合とで 、互いに大きさは異なる。即ち
ζij≠ζjjであり、さらに、上記式(0001)で述べた様に、Vij≠Vjiである。従って、bij≠bjiである。
またg(θ)はボンド角因子であり、
〔数34〕
g(θ)=1−(c2/d2)−c2/(d2+cosθ2) …(0006)
である。ここで、θは図4の様に取るものとする。θを求めるにあたり、実際の直交座標を用いて表現してみる。即ち、
〔数35〕
rij=√{(xj−xi)2+(yj−yi)2+(zj−zi)2 } …(0007)
であり、rikも同様の手続きで求められる。そうすると、内積をPijとすると、
〔数36〕
Pijk=(xj−xi)(xk−xi)+(yj−yi)(yk−yi)
+(zj−zi)(zk−zi) …(0008)
である。これらを用いて、
cosθijk=Pijk/(rij・rik) …(0009)
となる。ここで、本発明者等が用いた上記各式における定数を示す。即ち、
〔数37〕
R=3.0Å,D=0.2Å,A=3264.7ev,B=95.373ev,C=4.8381,
λ1=3.2394Å,λ2=1.3258Å,λ3=λ2 ,
β=0.33675,n=22.956,d=2.0417 …(0010)
である。
これらの準備を経ていよいよ本発明者の巡回部分を詳しく説明する。ここからいよいよ力を計算して行くわけである。ポテンシャルの(0002)式を位置の座標で微分すると力になる。即ち、
がそれぞれ、粒子i、jに働く力のベクトルのx成分である。しかし、実際にそれを求めるのは、そう簡単ではない。本発明者等の数学理論に立ち戻った工夫と発明が以下に有るわけである。それが、高階巡回偏微分と言う分けである。
(0011)式及び(0012)式で述べた式を、具体的に偏微分すれば良い。
しかし、そんなに単純では無いことが直ぐ分かる。
また、j及びkに関する偏微分方程式の変形は、同様に以下の様になる。また 、上記は主にxについて記述したが、yやzについても行う必要がある。ここでは、特に、上記との対応が分かるように、空白部分は、空白のままにして置いた。このような複雑な変形は、本発明者等が、調べた限りでは見あたらなかった。
この様な複雑な変形を今まで克服しなっかたため、逆に、単純なテーブル参照方式に逃げたりしていたわけである。この様な発明を克服したため、本発明者等は、分子動力学と汎用流体シミュレータとを合体出来たわけである。
これらの式の各項は次のように変形できる。
〔数41〕
∂Vij/∂rij
=∂fc(rij)/∂rij・[Aexp(−λ1/rij)−bij Bexp(−λ2/rij)]+
fc(rij)[−λ1 Aexp(−λ1/rij)+λ2bijBexp(-λ2/rij)]
=Aexp(−λ1/rij)・[∂fc(rij)/∂rij−λ2fc(rij)]−
Bexp(−λ2/rij)・[∂fc(rij)/∂rij−λ2fc(rij)]bij
…(0016)
∂fc(rij)/∂rij=−π/4D・cos[π/2・(r−R)/D]
(R−D<r<R+D)
∂fc(rij)/∂rij=0 (r≧R+D)
…(0017)
∂Vij/∂ζij=(∂Vij/∂bij)・(∂bij/∂ζij)
=−Bfc(rij)exp-(λ2rij)(−1/2n)(1+βnζij n)−1/2n-1βnζij n-1
=Bfc(rij)exp-(λ2rij)bij(βζij)n/[2{1+(βζij)n}ζij]
…(0018)
∂ζij/∂rij
=3λ3 3Σfc(rik)g(θijk)exp(λ3 3(rij−rik)3)(rij−rik)2
…(0019)
〔数42〕
∂ζij/∂rik=dfc(rij)/drij・g(θijk)exp(λ3 3(rij−rik)3)
−3λ3 3fc(rik)fc(rij)(rij−rik)2
…(0020)
∂ζij/∂cosθijk
=fc(rik)exp(λ3 3(rij−rik)3)dg(θijk/dcosθijk
fc(rik)exp(λ3 3(rij−rik)3)[2c2cosθijk/{d2cos2θijk}2]
…(0021)
∂rij/∂xi=(xi−xj)/rij=∂rij/∂xi
…(0022)
∂rik/∂xi=(xi−xk)/rik=∂rij/∂xk
…(0023)
〔数43〕
∂cosθijk/∂xi=1/(rijrik)・∂Pijk/∂xi +
Pijk{1/(rik)∂/∂xi[1/(rik)]+
1/(rij)∂/∂xi[1/(rik)]}
=1/(rijrik)・(xi−xk+xi−xj)−
Pijk{(xixj)/(rikrij 3)+(xi−xk)/(rikrij 3)
=1/(rik)・{(xi−xj)/(rij)−
(xi−xk)/(rik)cosθij}
…(0024)
∂cosθijk/∂xi=−1/(rij)・{(xi−xk)/(rik)
−(xi−xj)/(rik)cosθij}
…(0025)
∂cosθijk/∂xk=−1/(rik)・{(xi−xj)/(rij)
−(xi−xk)/(rik)cosθij}
…(0026)
L(ri∂ri/∂t,V,∂V/∂t)
=1/2Σm(∂ri/∂t)2−U(ri)+1/2M(∂V/∂t)2−PEV
(∂L(qj,q′j)/∂qj)−d(∂L/∂q′j)/dt=0
…(0027)
図7及び図8に、分子動力学シミュレータIIと、上記プロセス変動算出ユニットIとの連携を説明する。分子動力学シミュレータIIには、第一の入力aとしては、図8に示したファイル(22)、(32)、(42)、(52)、(62)、(72)等から転送されるものがある。内容は、境界条件や、領域の大きさ、応力、歪み、温度などである。また必要に応じて第二のものとして補助的な入力bがあり、分子動力学の計算時間やアンサンブルの指定等がある。
他方、第一の出力としては、位置ベクトルr(t)が出てくる。そして、図8に示したファイル(23)、(33)、(43)、(53)、(63)、(73)、(83)等へ転送される。叉、膜特性算出プログラムIIIにより、光学的・電気的・結晶学的特性等の膜特性を算出し、図8のファイル(3)へ転送される。
また、図9に、汎用2・3次元プロセスシミュレータ(図8の(Ia))の入力シーケンスの一例を示す。また、図10に、同じく汎用2・3次元プロセスシミュレータ(図8の(10))のもうひとつの入力シーケンスを示す。図9に示した内容は、例えば、それぞれ用いるクラスタツールをも指定するとともに、例えばその中にある酸化工程の詳細には、ガス量や、昇温シーケンス、さらには、バルブシーケンス等がある。
他方、図10の入力シーケンスは、図8で示した未知材料を含む場合である。
この図10の場合は、例えば、強誘電体膜であり、膜をスパッタする工程を含んでいる。このとき、強誘電体膜の詳細データがないので、分子動力学シミュレータIIに、一旦、飛び、そこで計算を行われる。
又、本発明者らは、経験的ポテンシャルを与えた分子動力学シミュレーションを行い、(1)各プロセス要素の揺らぎや各要素のずれの重畳現象も診断し、(2)所望のシーケンス通り実行されているかのチェックや、(3)変動因子の発見も可能 とするシミュレーションシステムを開発した。
具体的には、単結晶強誘電体膜とSi基板との位置関係や、金属電極との位置関連を提案するとともに、さらに、単結晶強誘電体膜として、実際上どの程度の品質のものを用意すれば良いかについても、初めて指針を示すことが出来た。かかる単結晶強誘電体膜構造を厳密に独自の計算手法により、計算機を用いて再現するとともに、特性評価関数を付与し、これに基ずき、印加条件を考えて基板単結晶との位置関係をどの様にすれば良いかと同時にその欠陥密度の許容範囲の指針を示すものである。
すなわち、本発明は、複数の工程からなる半導体装置の製造方法において、前記複数の工程の少なくとも一つにおける製造装置内に形成された膜、半導体基板上に形成された膜、半導体基板表面、及び半導体基板内部の少なくとも一つに関する複数の特性の少なくとも一つを測定して実観測データを得る工程と、経験的ポテンシャルを与えた分子動力学シミュレータにより、前記複数の工程の少なくとも一つにおける、複数の製造工程因子の少なくとも一つの設定値、及びこの設定ちの揺らぎを入力値とし、この入力値を基に変分計算を行い、前記複数の特性の少なくとも一つの予測データを得る工程と、前記予測データと実観測データとを、逐次、実時間で比較検定する工程と、前記比較検定により、前記複数の製造工程因子の少なくとも一つの設定値と、前記実観測データから推測される前記複数の製造工程因子の少なくとも一つとの間に有為差が認められた場合、前記製造工程因子を逐次実時間で修正処理する工程とを具備することを特徴とする半導体装置の製造方法を提供する。
本発明は、上述の半導体装置の製造方法において使用する経験的ポテンシャルを与えた分子動力学シミュレータは、前記複数の工程において形成または処理される半導体材料、絶縁材料、または導電性材料を構成する原子の運動が平行に達した所で、時刻を0に設定し、その後、時刻tまでをサンプリング時間とし、各原子の電荷及び位置ベクトルとの時間相関を求め、所定時間における光学的スペクトルを求める機能を有し、前記経験的ポテンシャルを与えた分子動力学シミュレータから導かれる予測データは、個々の原子の時時刻々の運動から導かれる振動挙動、その固有振動数、及び、固有振動数の分布情報を含む半導体装置の製造方法を提供する。
本発明は、上述の経験的ポテンシャルを与えた分子動力学シミュレータにおいて、原子の運動を計算する際に、原子の位置ベクトルと各元素の電気陰性度、もしくは各元素のイオン化ポテンシャルと電子親和力とを入力し、これらのデータから、個々原子の振動及び変位による電荷の変化を計算し、さらにこの電荷の変化を加味して前記ポテンシャルを計算する半導体装置の製造方法を提供する。
本発明者らは新しい計算化学理論により、種々元素を含む系における様々な現象を高速に且つ高精度に計算予測する方法を発明した。これによって、半導体製造プロセスに使用される様々な材料の特性が高精度に計算できるようになった。
従来、種々元素が含まれている結晶や分子の構造及び特性を予測する手段としては、abinitio分子動力学法や第一原理分子軌道法が使用されてきた。
abinitio分子動力学法は主に結晶の構造や電子状態の調査に使用され、分子軌道法は通常の分子や、結晶の一部を切り出したクラスタを対象として、その構造や電子状態の計算を行うのに使用されてきた。ところが、これらの方法では、系の各電子に対してシュレテ゛ィンガー方程式を解くため、膨大な計算時間が 必要であった。そのため、経験的ポテンシャルを与えた分子動力学法のように、数10万ステップもの計算回数を必要とする手法の中で、これら第一原理手法を取り込む事は現実的に不可能であった。一方、分子動力学の分野では、1991年に、A.K.RappeとW.A.Goddardが"Charge Equillibration for Molecular Dynamics Simulations"という論文の中で、多原子分子に対する新しい電荷計算手法を 提案している。本発明者らは、A.R.Rappeらによって提案されている多原子分子 における電荷計算手法をさらに大きく拡張し、種々の元素や結晶系に適用できるように大幅に変更し、半導体装置の材料として一般に用いられるSiやSiO2はもちろん、PbTiO3やPZTといった強誘電体や、高誘電体にも適用できる材料設計シミ ュレーション手法を構築した。まずはじめに、Rappeらの手法について以下に述 べる。
Rappeらは、まず、孤立原子のエネルギーEをイオン化ポテンシャルと電子親和力の関数として考えた。すなわち、ある原子Aにおいて、電気的に中性な状態をA0で示す。テイラー展開の2次の項まで考えると、
原子が電気的に中性な状態のエネルギーは、(3031)式のQAに0を代入して得 られる。
E
A(0)=E
A0 …(3041)
原子が+1価イオン化した時のエネルギーは、(3031)のQAに1を代入して得ら れ、-1価にイオン化した時のエネルギーは、-1を代入して得られる。
(3042)と(3041)の差はイオン化ポテンシャルIPになり、(3043)と( 3041)の差は電子親和力EAになる。
(3045)式のEAに−が付いているのは、電子を取るのを正の仕事と考えると、電子を加える事は負の仕事をした事になるためである。(3044)と(3045)の差をとると、電気陰性度χA0となる。
また、(3044)と(3045)の和をとると、次式が成り立つ。
ここで、(3052)式の物理的意味を考える。水素原子の場合のイオン化時のクーロン相互作用の模式図を図12乃至図14に示す。図から分かる様に、IPとEAの差は電子間クーロン相互作用bになっている。従って、φA軌道中の2電子間 のクーロン斥力をJAA0と表すと、
IP−EA≡JAA 0 …(3053)
JAA0をidempotentialと呼ぶことにする。電子を加える事により、軌道形状が 変化するので、Hartree-Fock等の第一原理計算で計算されたJAAと(3053)式 の値は多少異なってくるが、ここでは(3053)式で近似する。
Idempotential JAA0と原子直径RA0の関係を考えてみる。図15に示した簡単 な原子構造モデルを考えると、電子間のクーロン相互作用JAA0は以下の様になる。
ここで、 JAA0はeV単位、RはÅ単位である。(3053)式の妥当性を確認するため、J値とRとの関係を調べ、表1に纏めた。表1に示した様に、Jは原子サイズに ほぼ逆比例し、(3054)式から求められた原子の直径は等極結合距離とおおよそ一致していることが分かった。
(3051)(3053)式を用いると、孤立原子のエネルギー(3031)は、次式で書ける。
多原子分子の場合には、原子間の静電エネルギーも考慮する必要がある。静電エネルギーは次式で表せる。
ここで、JABは原子A,Bの中心に1電子電荷がある時のクーロン相互作用である。JABはAB間の距離に依存する。(3055)式を用いると、全静電エネルギーは(3 057)式で表せる。
原子スケールの化学ポテンシャルχI(Q1,・・・・,QN)は次式で表せる。
電荷平衡の条件(3059)から、(N-1)ヶの連立方程式が得られる。
χ
1=χ
2=………=χ
N …(3059)
すなわち、
となる。これに、全電荷の条件(3062)式を含めて、N個の連立方程式を解く 。
(3061)(3062)をマトリックスで表示すると、以下の式で表せる。
(3063)式を解くに当たって、各イオンが取り得る電荷には、上限下限があるので、電荷がその範囲外ならば、その電荷を境界値に固定する。(3061)式を境界値に固定した電荷と、そうでない電荷とに分離し、(3063)式を修正してマトリックスを解き、収束させる。すると、電荷のばらつきが計算できる。
この方法の中で問題なのはJABの決定である。原子AとBとの距離RABが大きいときには
でよいが、電子雲が重なる距離においては斥力を考える必要がある。ここでは、1つのスレーター軌道の項で電子密度を近似してしまう。すなわち、外側の原子 価軌道がns,np,あるいはndである原子において、何れの場合にも次の型の規格化されたnsスレーター軌道からなるとする。
(3065)式を用いて原子の平均サイズを求めると、
となる。下記の関係によって、原子価軌道指数ζAを選ぶ。
λは調整用のパラメタで、(3066)式で与えられた平均的な原子サイズと結晶の共有結合半径RAとの差を計算するために含まれている。
スケーリングパラメタλは、図16に示したようなアルカリ金属ハライド分子を用いて設定する。アルカリ金属ハライド分子をMXとし、Mをアルカリ金属( Na、K、Rb、Cs)、Xをハライド(Cl、Br、I)とする。それぞれの電 荷をQM、QXと表すと、(3062)式から、
Q
total=Q
M+Q
X=0 ∴Q
X=−Q
M …(3068)
(3061)式から
(3069)式に(3068)式を代入すると、
〔数62〕
(J
MM−J
MX−J
XM+J
XX)Q
M=χ
0 x−χ
0 M …(3071)ところで明らかに
J
MX=J
XM …(3072)
であるから、
(3073)式の変数はλだけである。計算されたQから(3074)式を用いて双 極子モーメントμMXを求め、実験値と比較してλを決定する。
μMX=4.80324QMRMX …(3074)
RMXは結合距離の実験値である。(3074)式はQがe単位、RがÅ、μがデバイの単位になっている。双極子モーメントの実験値との比較から、λの最適値として0.4913が得られている。このλを用いると(3067)式より原子価軌道指数ζが図17の表の値に求まる。
水素に関しては、Rappeらは、他の元素と異なる取り扱いをしている。
水素の電気陰性度をMullikenとPauling他の実験上の値と比較すると一致して いない。Paulingスケールでは、水素は炭素より電気陽性で、ホ゛ロンよりわずかに電気 陰性であるのに対し、Mullikenスケールでは、水素は炭素や窒素より電気陰性である。この電気陰性度の問題は結合している水素の軌道は自由な水素イオンのように広がる事ができないので、実行的な電子親和力EAが原子の時の値よりも小さくなるためである。
従って、EAHを変数にして、水素のχ0HとJ0HHとを再定義する。
実験との比較から、最小二乗法によって、
χ0 H=4.528eV
J0 HH(0)=13.8904eV …(3076)
実験でなく、HF波動関数の静電ポテンシャルから計算された電荷との比較から、 χ0 H=4.7174eV
J0 HH(0)=13.4725eV …(3077)
を使用しても良いとしている。
また、電荷に関する境界値を設定する。
χA 0,JAA 0は、明らかに、原子価殻が空っぽ、あるいは満杯に対応する範囲 外では妥当でない。例えば、Li,C,Oの場合には、下図に示したような範囲に電荷を制限する。この範囲外ではEA(Q)=∞をとる。
マトリックス(3063)を解いて求めた電荷を、図18乃至図20に示した不等式の範囲内にあるか否か判定し、その範囲外である場合には、その電荷を境界値に固定する。
以上がRappeらの手法のあらましである。Rappeらは、前述の手法を有機系の分子に適用している。本発明者らは、Rappeらの提案した原子スケールの化学ポテンシャ ルという概念だけを採用し、結晶系における個々原子の電荷計算手法を新たに開発した。以下に本発明者らによる電荷計算手法を述べる。
(3057)式右辺第二項は、すべての結合の静電エネルギーの和を意味する。本発明者らは、分子から結晶に拡張し、周期境界条件を取り入れようと試みた。すると、静電エネルギーを計算する際、計算に用いている数百〜数千の原子だけでなく、その外側にある無数の仮想原子からうける静電エネルギーも考慮する必要があることを発見した。すなわち、静電エネルギーは減衰が小さく、遠方の原子からの影響も蒸し出来ないことに気がついた。そこで、計算に用いているセルの外側に仮想的に無数のセルが存在すると仮定して、仮想セルからの静電エネルギーも考慮するようにした。
この場合、(3057)式仮想セル中の原子から受ける静電エネルギーJABは 、計算セルの大きさを考慮すると、RABが通常10オングストローム以上の大きな値 となるので、クーロンエネルギーだけと考がえて良い。従って、JAB=1/RABと一 律に扱う。仮想セルをνで表すと、(3051)式は次の様に変更される。
ここで、ν=0は計算セル自体を意味し、ν≠0は仮想セルを意味する。(3078)式 第3項にエワルト゛法を適用する。エワルト゛法を用いると、第3項は次の様に求めることが出来る。
Gは逆格子ベクトル、Ωは計算セルの体積、κはエワルド計算の収束を調節する パラメタである。更に、(3082)式右辺第一項はν=0且つA=Bの項は含まず、 第二項はG=0の項は含まない。これらの式から、原子スケールの化学ポテンシャルχIを求めると、
となる。電荷平衡の条件から、N-1個の方程式が得られる。
χ
1=χ
2=χ
3=……=χ
N …(3092)
(3091)式を代入して、
(3093)式で表せるN個の連立方程式を解くことになる。
(3093)式からCQ=-Dのマトリックスを作ると、次式が得られる。
本発明者らは、SiO2の様な結晶や非晶質を扱う場合には、周期境界条件を考慮しながら(3094)式を解けば良い事を発見した。(3094)式は電荷Qに対し て非線形となっているため、計算プログラムにはQを収束させるためのルーフ゜を持 ってくる。
さらに、本発明者らは、実計算にあたってのクーロン積分のモデル化にも注意を払った。JABは、ζAとζBに対するスレーター関数の2原子クーロン積分として求め た。クーロン積分JAB。はRoothernの方法を用いて計算した。Roothernらは、ク ーロン積分JABには、 という標識を使用する。 (3101)式で定義されたハ゜ラメタを用いて、クーロン積分は(3102)式で表せる。
ここで
〔数71〕
(α,α′,β,β′)
=(1+τ
a)α(1−τ
a)α′(1+τ
b)β(1−τ
b)β′
…(3103)である。さらに、
ここで、上式の単位系は原子単位系であり、距離はbohr、エネルギーはHartreeを使用している。
Roothernらは1s,2s,2pに対するクーロン積分の解析式を与えているが、本発明者らは、3s以降の原子軌道に対するクーロン積分まで計算し、半導体装置で一般的に用いられるSi及びSiO2が扱えるように、最適な軌道の選択を行った。
JABの値は電荷計算に大きく影響する。そのためJを算出する際に用いるスレーター軌道や原子価軌道指数ζの設定は重要である。そこで本発明者らは、Si-OやSi-Siを扱う場合のJの値について調べた。
Rappeらは「分極電荷がns軌道にあると仮定」しているが、種々元素に対して どの軌道を考慮するのが最適であるかという点は不明である。 Siの主量子数nが3であることから、Siの電荷が3sにあるとモデル化した場合と、2sにあると仮定 した場合の計算を行い、比較した。このとき、n=3におけるJの解析式はRoothernの方法に従い、本発明者らが、別途用意した。例えば、Si-Siに対しては、次式 を用いた。
図21はSi-Si間のクーロン相互作用Jの計算結果である。この場合には、2sと3sの何れを使用してもJの値に殆ど変化が無いことが分かる。ところがSi-O間の クーロン相互作用Jは2sと3sを使用した場合で大きく異なっている。図22にOの電荷が2sにあり、Siが2sあるいは3sにあるとモデル化した場合のJABの計算結果 を示した。Siの電荷を3sだけにおくと、3s軌道が非常に広がっている影響を受けて、Jの値がr=1.6Å近辺で極小値を持ってしまう。Si-Oの結合距離は丁度この当たりにあるので、Jが非常に小さく見積もられ、結果、SiとOの何れの電荷も極めて過小に評価されてしまう事が分かった。Si-Hの場合も同様であった。従って、Siの電荷のモデルとしては、敢えて3s軌道で仮定すると、Jの見積もりが全く精 度を欠いてしまう。このことから、スレーター軌道を用いるのは、電荷の広がりを考慮するためだけとし、モデルとしては、1sあるいは2sのスレーター軌道を用いるのが最適である事が分かった。これは、Si以降の原子番号の元素に対しても同様であり、3s以降の広がった軌道をモデルとすると、原子間距離が近い場合のクーロン力の計算精度が非常に悪くなることが分かった。
水素に対する補正にも、本発明者らは、プログラム上の注意を払った。 (3031)式から、
水素補正した時の原子スケールの化学ポテンシャルχHlを求める。まず、(30 57)式に対応する多原子分子のエネルギーの式は、全原子数N個の内、水素M個 とすると、
ここで、(3114)式の右辺第一項、第二項、第五項のΣ′の記号は、それぞれ の和を求めるに当たってHは除いている事を意味している。原子スケールの化学 ポテンシャルχHlは、
(3121)式右辺第四項は次の様に解きほぐす事が出来る。
(3122)式より、(27)式の第四項と第五項との和は、Hl以外の全原子との間の原子間クーロン相互作用を意味する事が分かる。従って、化学ポテンシャルは、
(3123)式を用い、水素補正時のマトリックス(3063)に変更を加える事になる。今、χ1が水素以外の原子であるとすると、(3059)式に(3123)式 を代入し、
マトリックス(3063)と比較すると、水素のJii0を
に変更する事になる。
次に原子間クーロンポテンシャルJAHの水素補正を考える。水素の場合のスレーター軌道関数は、(3075)式より、
水素が含まれる場合には、電荷を求めた後に(3126)式を用いてJAHを修正 し、再度電荷を求める。そして電荷の値が収束するまで、これを反復する。この際、I=1の行が水素の場合、水素でない行と入れ替えてマトリックスを作成する ように、本発明者らはアルゴリズムを作成した。
電荷境界条件の扱いにも、本発明者らはマトリックスの具体的な変更方法について注意を払った。
(3061)式を、境界値に固定した電荷QBと、そうでない電荷QCに分けて考える。
(3062)式を(3131)式に変更すると、マトリックスの要素は、のとき、
上式の様にマトリックスを変更して、電荷を境界条件の範囲内に収める。この時、I=1の元素が境界値に設定される場合、マトリックスの行を入れ替え、境界 値設定されていない元素に対する行がI=1になるように行の選択(ピボット選択) をするようにした。
本発明者らは、結晶系及び非晶質系に新たに開発した上記電荷計算法を用い、これを分子動力学法に追加して、個々原子の動きを逐次計算し、双極子モーメントを時々刻々計算して記録し、これをフーリエ変換してIRスペクトルの算出を行った。図23は、以上の手続きをまとめたものである。
一方、16G以降のDRAMに使用されるゲート酸化膜やNAND EEPROMのトンネル酸化膜においては、高信頼化が従来にも増して重要な課題となっている。特に、EEPROMにおいては、高電界下でトンネル電流を流しながら絶縁性を維持させるという過酷な条件下で酸化膜を使用するため、絶縁材料としての限界を原子・電子レベルで探索することが重要となっている。
このような状況から、Si/SiO2界面構造やSiO2ネットワーク構造の原子レベルでの探索実験や、SiO2のトラップ機構の実験、及び理論計算による 予測等が種々の研究機関で進められている。しかし、本発明のように、ネットワーク構造とIRフルスペクトルとの関連を原子レベルで関連づけ、利用する手法は、従来にはなかった。僅かに、IRスペクトルのSi−Oの非対称伸縮ピークと酸化界面の応力とを関連づける実験結果が報告さているだけであった。
また、実験結果の解釈としても、「中心力モデル」を用い、現実のネットワークを非常に簡略化したSi−O−Si分子構造の振動解析からの議論しかなされておらず、非晶質特有のSi−O結合長の分布、Si−O−Si、O−Si−O結合角度の分布等を全て考慮した議論は、従来全く行われていなかった。本発明においては、非晶質特有のこれら様々な分布を考慮することにより、中心力モデルとは異なった視点で、電気的特性と構造との相関が得られたのである。
ここで、本発明者らが克服した点について付言する。本発明者らは、古典的分子動力学法を使用するにあたり、そのポテンシャルを吟味した上、さらに、改良を加えた。本発明者らは、最初、常行氏らによって開発されたTTAMポテンシャルの使用を試みたが、Si−O非対称伸縮振動が非常に低波数側に表されるという欠点を見い出した。そこで、van Beest等によって開発されたBSKポテンシャルを使用した。彼らのポテンシャルは、TTAMポテンシャルよりも高波数側にピークが表れたが、尚、実測値よりも低い波数に伸縮ピークが表された。
本発明者らは、このピーク値を実測に近づけるため、これら既存ポテンシャルで使用されている、電荷一定という条件を取り除き、Si−O−Si角等が大変位した場合の電荷の分布を、古典的分子動力学に取り込んだ。この電荷分布の計算は、理想的には分子軌道法を用いて算出できるが、現実問題としては、数百原子もの構造を一度に取り扱って電荷分布計算を行うことは現状では不可能である。スーパーコンピュータを用いても現実的には不可能である。
そこで、本発明者らは、A.R.RappeとW.A.Goddard III によって提案されている手法を用いて電荷分布計算を行った。引用論文は“Charge Equillibration for Molecular Dynamics Simulation”,J.Phys.Chem.,95,1991,p.3358−3363である。本発明者らは、彼らの論文を参考に数式を展開し、論文中の誤りを訂正の上、必要なパラメータ算出を独自に行って、彼らの手法をSiO2系に適用した。Rappeらの手法の概略について以下に 述べる。
原子iの個別の電荷によるエネルギーEi (Q)を、テーラー展開の2次まで考えて、下記式(1651)で近似する。
ここで、Xi 0 は電気陰性度であり、下記式(1652)で表される。
Xi 0=1/2(IP+EA) …(1652)
Jiiは、原子Iの軌道上の電子間のクーロン斥力であり、下記式(1661)で定義される。
Jii=IP−EA …(1661)
ここでIPはイオン化ポテンシャル、EAは電気親和力を示す。N個の原子で構成される分子の全エネルギーE(Q1 ,…,QN )は、個別の原子のエネルギーと原子間の相互作用静電エネルギーJ
ijQi Qj の和である。ただし、J
ijは原子
ijの中心にある単位電荷のクーロン相互作用である。従って、下記式(1662)が成立する。
原子スケールの化学ポテンシャルは、下記式(1663)で与えられる。
この式(1663)を、電荷平衡の条件、X1 =…=XN に代入すると、下記式(1671)が得られる。
さらに、下記式(1673)に示す全電荷の条件を考慮する。
J
ijは既知であって、およそ下記式(1674)である。
ここでRはI=jのときは原子の大きさ、I≠jのときは原子間距離である。従って、行列Cとベクトルdを導入すると、Qに関する下記式(1681)に示す連立方程式が成立する。すなわち、
と定義すると、
CQ=−d …(1682)
となる。論文ではC1i=Qi 、Cd =dとあったが、本発明者らは、これらを上述の様に修正して使用した。さらに、Si−O−H等の原子の大きさ等、必要なパラメータを独自に算出し、連立方程式を解いて電荷分布を求めた。これによって、ピーク波数が格段に実測と一致するようになった。
分子動力学シミュレータIIで算出された「戸籍簿」、つまり単位時間毎の結合長(位置ベクトルr(t))と結合角を示す膨大なデータの一例を、図24に示す。このデータは、膜特性算出ユニットIIIに入力され、ここで基本膜特性(光学的、電気的、結晶学的特性)が算出される。
又、プロセス変動算出ユニットIについては、その2次元汎用システムの一例が、特開平3−164904号に記載されている。
3次元のプロセスシミュレータは、今回本件発明者が構築に成功した。以下これを、図8をもちいながら説明する。まず技術とシミュレータの構造の概観をする。本願は、原子レベルに主体を置いているが、他のプロセス変動算出ユニットIとの新しい合体方法をも提案しており、より効果的かつ正確で、従来にない新しい機能をも提供できることをも示している。合体する相手のプロセス変動算出ユニット(10)は、例えば3次元プロセスシミュレータIcであり、いわゆる流体モデルに属するものである。とくに、本発明者らが新しく構築した3次元プロセスシミュレータIcの内容を以下に記す。要点は2つあり,まず、本発明者らが構築した「3次元プロセスシミュレータ」は、原子レベルシミュレータとの接続が可能であること。次に、この「3次元プロセスシミュレータ」の部分には、本発明者等が新しく克服した化学モデルや、数学式を駆使している。これらにより、初めて、原子レベルとの合体が可能になったのである。原子レベルシミュレータから出力されるものとして例えば、粘弾性係数、応力定数、拡散定数、Cij、等があり、これらが、流体的汎用3次元プロセスシミュレータに転送されるものである。図8では、(23)、(33)、(43)、(53)、(63)、(73)、(83)である。例えば、酸化剤等の拡散定数は、(23)に転送されるわけである。転送手続きの方法等は後述する。
汎用3次元プロセスシミュレータ(10)側で、入力データに従って、数ミクロンレベルの広領域の工程誘起応力の計算を、(2)や(3)で、行い、この結果を、ファイル(22)やファイル(32)を介して、一旦分子動力学シミュレータIIに転送し、あるいは、再び汎用シミュレータ(10)に戻し、これを繰り返すことにより、より正確でより多くの機能を利用することもできる。計算手続きとしては、例えばある時刻で、プロセス変動算出ユニットで求めたある程度広領域の工程誘起応力計算の結果を、一旦分子動力学シミュレータIIに転送し、ここで、物性定数の修正量を計算しながら、改めて、原子の動きをみ、再び、修正された物性定数を汎用シミュレータに送り込み、汎用シミュレータで広領域の工程誘起応力等の結果を、進める。これを繰り返すことにより、もし、十分大きな応力が局所的に生成されれば、分子動力学シミュレータII側で原子の動きをみ、塑性変形の可否を見ることもできる。
勿論上記の汎用プロセスシミュレータ部分での不純物再分布計算や、応力計算の時には、図8に記した、点欠陥挙動計算部分(9)が呼ばれている。点欠陥としてはSi基板中の空格子と格子間Siを扱っている。不純物の拡散にはこれら点欠陥の挙動を正しく取り入れないと、不純物の再分布は正確には計算されない。イオン注入計算部分(4)では、もちろん、注入損傷による多量の空格子と格子間Siが生成される。これらは、熱工程中に再分布するが、その時、不純物とも相互作用を持ち、不純物分布の計算に反映させなければならない。またシリサイド・形状・応力計算部分(7)でも、勿論、Si基板中の空格子と格子間Siの挙動に十分配慮しなければならない。
ところで、図8のその他(8)についても、少し触れる。LSI素子試作工程で、酸化・拡散・イオン注入等は、既に多くのデータを基に、このプロセス変動算出ユニット(10)には、各必要物性定数が格納されている。例えば、Siの酸化に関しては、後に掲げる式(1130)や(1131)のPB1やPB2等がそれである。また、SiO2膜中の酸化剤の拡散定数などは、後に掲げる式(1531)のD’などがそれである。一般には、この様に、概ねデータがあるものである。
しかし、LSI開発の中では、必ずしもデータ等が揃っているとは限らない。
たとえば、今、強誘電体膜を用いたメモリーセルの試作開発を考える。CDV或いはスパッタ工程を用いて、PbTiO3を形成する場合を考える。現状の汎用 プロセスシミュレータでは、多くはまだPbTiO3膜中での種々の不純物の拡 散係数も、さほど分かっていないし、また、その応力定数もさほど分かっていない。また、界面での不純物の偏析定数もさほど分かっていない。この様な場合、まず(8)の工程に入り、すぐさま、(84)に示す様に、分子動力学シミュレータIIに一旦行き、そこで、必要な元素のポテンシャル等を算出し、それを基に、分子動力学部分が動き、そして、物性定数を算出する。このあとは、既に述べてきた部分と非常に良く似ているが、(83)からこれらの物性定数がファイルに格納される。そして、プロセス変動算出ユニットと同じようにして、計算部分が進められる。
ここで、図9を参照して、具体的な処理に即した説明を行う。ここでは、酸化工程を取り上げる。即ち、ウエハーの投入を行い、これを基板処理工程に移し、しかる後に、酸化工程に入る。
この酸化工程では、例えば、高純度酸素ガスの量を例えば、Poxとし、また、 これの分散をσtoxとして構成する。叉、昇温シーケンスはTempとし、またこれの分散をσtemp、バルプシーケンスをtox、σtox等の様に構成する。特に、バルプシーケンスは、ここでは、例えばパルプを開けている時問をさしているわけであるが、詳細は実施例で述べたようなものである。
たとえば、図8では、(2)の酸化・形状・応力計算部分をアクセスし、ここに分圧Poxや、温度Temp等が入力される。そして、酸化成長計算や、同時に応力計算も進められる。この時、(21)を用いて、分散・変動計算も行う。ここではσpoxやσte即が入力して、中心値からの変分を効率よく計算するものである 。
これらの計算結果は、例えば、酸化膜の成長量や、酸化膜とSiの界面の形状や、さらには、不純物量も計算される。さらに、本発明者によれば、応力も求められるわけである。また分散・変動からも同じように Pox=Pox+σpoxや、Temp=Temp十σtemp等の部分についても、酸化膜の成長量や、酸化膜とSiの界 面の形状や、さらには、不純物量が計算される。そして、これらが、一旦メモリまたはファイル(22)にたくわえられ、そして、分子動力学部分に送られる。
分子動力学部分では信号aとしてこれらが入力され、位置ベクトルr(t)が算出される。
図10は、図9と同様の行程図を示しているが、ここでは強誘電体膜の形成に関するものである。
ところで、ここで本発明での診断部分に関する説明も行っておく。図11に、本実施例による2進木の一部を示す。ここでは酸化工程に着目した2進木を示している。図中にも記したように、ここではまず、in-situモ二夕から送付される SiO2膜厚と、計算から予測されるSiO2膜厚を比較する。そして、膜厚が許容範囲に収まっている時を取り上げ議論している。そうすると、次の手続きとして、inlituモニタから送付される1Rスベクトルカープと、計算から予測される1Rカーブとを比較する。そして、図に示した様に、まず、スペクトルを比較して、少しずれていた場合で、しかも、スペクトルが高波数側にずれた場合を取り上げる。そうすると、次の手続きとして、SiO2中に残留応力が大きく誘起し ていることが考えられる。それで、この場合、ゆっくり降温して、残留応力を軽減させることが出来る。それで、図9に記載したアニール部分を修正して、図8の(2)に入力される。他方、スペクトルを比較し、許容範囲に入っていた場合は、当初のプロセスシーケンスのまま続行される。
図8に関して、概ねデータの流れがこれで分かったと思う。では、もう少し、図8の出力に関して言及したい。即ち、実行結果は、一旦、ファイル(3)に蓄えられる。ここでは、各工程で求めた、工程に対応するスペクトルや、その他の光学的・電気的特性結果が格納される。これらは、あとで、詳細を述べるが、実測データとしてモニタをしている量と対比できる様になっている。
一方、本発明者らの1人は、以前、プロセス変動の変分計算を用いたプロセスシミュレータシステムを提案している。即ち、特開昭63−174331号および特開平3−164904号である。これらの提案では、流体的な取扱いにより、酸化/拡散、イオン注入、化学的堆積(CVD)、エッチング、リソグラフィ等の各工程を取り扱うことができる。例えば、流体モデルによる拡散工程では、不純物総量と電気的活性化不純物量とを区別しながら、電気効果を取り入れつつ、不純物の再分布を求めている。
特開昭63−174331号および特開平3−164904号に示すシミュレーションシステムの概略を、図25〜図49を参照して説明する。これらの提案では、プロセス変動については次のように扱っている。例えば、酸化、イオン注入、リソグラフィの工程を例にとると、平均の温度(T)とその許容バラつき(δT)、平均時間(t)とその許容バラつき(δt)、平均圧力(P)とその許容バラつき(δP)、更に平均マスク寸法(L)とその許容バラつき(δL)、平均ドーズ量(D)とその許容バラつき(δD)などが入力値である。
プロセスシミュレータには、図25に示すように、イオン注入、酸化/拡散、リソグラフィ等の各サブルーチンに分散変動(deviation)を計算する部分が付設されている。計算手順としては、これらの中心値について計算を進めながらその変分を計算して行く方法を採用している。
この変分計算について若干説明する。まず拡散工程における温度変分について説明すると、拡散工程では、拡散係数をD、不純物分布による内部電界をε、電気的活性化不純物数をNとして、Fickの第2法則に基づく方程式を解いて、不純物濃度分布Cを求めている。図26を用いて変分計算の手順を示す。ここで、図中の左側を正常プロセス(Normal process)とする。
非線形拡散方程式は、下記式(0071)により表わすことが出来る。
ここでk、k+1は時刻であり、時刻k+1のときの不純物濃度Ck+1 が未知数である。右辺のBは、下記式(1072)の第2項に由来するものである。
∂C/∂t→(Ck+1−Ck)/t …(1072)
ここでは陰解法を用いているので、上式に現れる係数行列はすべてDや平均温度T等を用いて表現されている。正常プロセスの場合は、所定の拡散時間tn+1 番目までを適当に分割し、t1 ,t2 ,…,tn+1 まで順次解を求めて行く。ここではそれぞれの時刻、t1 ,t2 ,…,tn+1 においてNewton法を用いてこれらを線形化して計算をすすめている。
次に、変動プロセス(Process deviated)について説明する。ここでは拡散温度に変動δTが発生した場合を示している。上式の変分をとり、時刻kにおける濃度の変化量δCK がわかっていると、時刻k+1の時の変分δCk+1 は、下記式(1081)により表される。
すなわち時刻k+1において正常プロセスの収束時の値D,ψ,N,C等を保存しておき、これらを用いてF′を作成し、またその温度Tに対する変化量を用意しておき、δCk+1 を求める。ここでF′を作成するための行列要素の一部を図26に示す。拡散係数の導出には、基本的には空孔モデルを用いている。この場合、温度変化によってフェルミレベルが変化し、これによって荷電空孔の濃度が変化し、拡散係数が変化する。従って、まずその基になる真性キャリア濃度ni に対してはその変分は、下記式(1082)に示すようになる。
〔数94〕
δni=ni(1.5/T+0.605/kT2)δT …(1082)
また、不純物の内部電界の温度による変化も、これらを用いて表現できる。更に、酸化性雰囲気における拡散係数の増大現象(Oxidation Enhanced Diffusion Effect)の変分に対しても考慮した。酸化工程では線形酸化速度定数(Linear rate constant)及び放物線酸化速度定数(paraboric rate constant)も変化するのでこれも考慮した。これらを用いて先のF′を作成する。
ここで、ある典型的な拡散工程を取り上げ、その1つの拡散時間における収束の様子を図27に示す。図27の横軸は反復回数を示し、縦軸は残差を示している。図27では、正常プロセスにおける収束状況が示されている。正常プロセスでは、3回のNewton loopで収束していることがわかる。この収束時のD,ψ,N等を用いて先に示した変分行列F′を作成し、温度変分を計算する。図27に示すように、約50回で収束し、全体の1/10の計算量で済んでいることがわかる。
また、図28、図29、図30は、変分の精度を示す。図28は、正常プロセスとして1000℃の時の不純物分布を2次元格子にわたって表現したものである。図29は、変動プロセスであって、δT=1℃の場合で図28をもとに計算したものである。他方、図30は、1001℃を正常プロセスとして打ち出したものである。
これらの図から、図30の値=図28の値+図29の値なる関係があることがわかる。また、改めて図30の値をプロットしてみると、図31に示すようになり、高温にズレると高濃度側では不純物濃度が下がる傾向にあり、低濃度側では増大する傾向にあることがわかる。
半導体素子を作成する場合、多くのプロセスがあり、各プロセスにおいては複数個のパラメータがある。例えば、CMOS工程では主なものだけでも、(i)N−well領域へのイオン注入、(ii)P−well領域へのイオン注入、(iii)LOCOS工程のためのSi3N4のパターニング、(iv)LOCOS工程、(v)n+ イオン注入、(iv)p+ イオン注入、(vii)アニール工程などからなる。プロセスの揺らぎを考える場合、これらの数多くのパラメータの変動とその度数(或は頻度)を取り扱う必要がある。
特開昭63−174331号および特開平3−164904号では、変分の度数と実行の組み合わせについて考えている。今、簡単のためにα,β,γの3つの工程を考えてみる。それぞれにおいてプロセスパラメータの中心値をC(Center)、正に変位した値をU(Upper)、そして負に変位した値をL(Lower)と名づけ、それらの度数を、特開昭63−174331号および特開平3−164904号から、図32に示す様に定義してみた。
工程αとβとを考えた場合、組み合わせは、UUやLLなどの計算実行を無視すれば、CU,UC,CC,LC,CLの5個の組み合わせが考えられる。さらにγのプロセスを引き続き経るときは、図中に示した7つの組み合わせが考えられる。しかし、特開昭63−174331号および特開平3−164904号では、これらを全て取り扱うのではなく、UとLの対称性を考慮し、Lの場合はUから計算するようにしている。
このシステムの中心部である推論エンジンの探索の木の一部を、図33に示す。ここでは、酸化・拡散工程について検討している。図中の“予”は、統計解析システムによる予測値を示し“実”は、実測値を示す。
このとき用いる許容値は、酸化・拡散・エッチング等全工程を経た値である。
このときε値は、2次元プロセス統計解析システムを実行して得られるものである。“実”が“予”±εの外に出た場合は、図33に示すように、酸化工程に異常があった旨記述する。この場合、温度又は時間がどれだけズレていたかを必要に応じて、統計解析システムのロードモジュールを実行させる。上記検定は、酸化・拡散工程を中心に説明しているが、リソグラフィ工程についても行う。リソグラフィ工程でのゆらぎはSiO2膜の形状の横ズレ、接合形状の横ズレなどに 相当する。
推論エンジンには、Common Lispを用いているが、これは数値処理のみならず文章(List)の処理、形状の処理に適していること、さらにシステムの保守が極めて簡単であるためである。本システムは、推論エンジンから実測値のデータを読むだけではなく、統計解析システムのロードモジュールをも実行している。
特開昭63−174331号および特開平3−164904号に示すシミュレーションシステムによると、各工程における正常反復計算の収束時の行列の係数を保存しておき、且つこれを用いて一次変分行列を作成し、前記個々のプロセスパラメータの変動、及びこれらの重畳効果をも考慮することにより、従来の約1/100の計算量で変動部分が処理できるようになった。また、たとえばCommon Lispによる推論エンジンと上記システムを連動させることにより、プロセスマージンさらにはプロセス異常の診断さらに最適化まで出来ることがわかった。
特開昭63−174331号および特開平3−164904号では、基本拡散方程式は、各不純物に対して下記式(1111)に示すように立てている。
ここで、x,zは2次元空間の水平方向と鉛直(下向き)方向の座標であり、tは時間を示し、Cは不純物濃度、Nは電気的活性化不純物量である。またDは濃度依存性と増速拡散効果を含む拡散係数であり、下記式(1112)により表わすことが出来る。
式中、Di は濃度に依存しない基本拡散係数であり、その温度依存性は、下記式(1113)で表現される。
D
i=BB・
exp(−BA/kT) …(1113)
BBとBAは不純物の種類による定数、kはボルツマン定数である。FEはni を真性キャリア濃度として、下記式(1121)で与えられる。
またnは1/2[(CD −CA )+(CD −CA )2 +4ni 2 ]である。さらにni は、
〔数98〕
ni=3.87×1016×T1.5×exp(−0.605/kT)
…(1122)
とした。
酸化膜成長速度は、放物性酸化速度定数R1と線形酸化速度定数R2とを用いてそれぞれ下記式(1123)、(1124)に示すように表現した。
R1=PB3・exp(PB1/kT)
…(1123)
R2=PB4・exp(PB2/kT)
…(1124)
上記放物性酸化速度定数をあらためてB、B/Aと記すと、水平方向の座標xの点で、酸化膜厚をZo 、酸化時間をt、τを定数として、下記式(1125)に示す関係がある。
Zo=A/2・[1+4B(t+τ)/A2−1]
…(1125)
これらを用意しておき上記各係数の温度δTに対する変分を求めてみる。
δni =ni ・(1.5/T+0.605/kT2 )δTとしてnの変分は、下記式(1131)に示すようになる。
成長係数の変分は、下記式(1133),(1134)により表わされる。
B=R1,A=R1/R2に書き改めると、下記式(1141),(1142)に示すようになる。
そして、t+τの変分は、下記式(1143)のようになる。
以上説明したシミュレーションシステムは、当時としては画期的なものであった。しかし、最近のクラスタ化ツールにおける工程制御やモニタ技術には、もはや使えなくなっている。
そこで本発明者等は、更に検討を重ね、本件発明を完成するに至った。ここでは、ペロブスカイト単結晶を用いた電界効果型MIS素子の強誘電体膜の最適設計を行い、これに従い、実際に素子を試作し評価し発明完成を確認した。
即ち、請求項11によって提供されるのは、半導体単結晶基板上に電界効果型MIS素子を形成するにあたり、該MIS素子の強誘電体膜を単結晶とし、しかも、該半導体基板と該単結晶強誘電体膜を含む系に於いて、対象特性を表現するabinitio分子動力学理論に基ずいた評価関数を導入し、これにより所望の主動作条件下での系の最適解に従った単結晶強誘電体膜および基板の設計配置とすることを特徴とする電界効果型MIS半導体装置の製造方法である。
叉、請求項12によって提供されるのは、半導体単結晶基板上に、電界効果型MIS素子を形成するにあたり、該MIS素子の強誘電体膜を単結晶とし、しかも、該半導体基板と該単結晶強誘電体膜を含む系に於いて、その系の対象特性を表現するabinitio分子動力学理論に基ずいた評価関数として系全体の自由エネルギを取り上げ、主動作条件下で自由エネルギが最小になるように前記単結晶強誘電体膜を配置することを特徴とする電界効果型MIS半導体装置の製造方法である。
更に、請求項13によって提供されるのは、Si半導体単結晶基板上に、電界効果型MIS素子を形成するにあたり、該MIS素子のゲート絶縁酸化膜を単結晶とし、しかも、該半導体基板と該単結晶強誘電体膜を含む系に於いて、その系の対象特性を表現するabinitio分子動力学理論に基ずいた評価関数として系全体の自由エネルギを取り上げ、この自由エネルギが、主動作条件下で、最小になるように前記単結晶強誘電体膜を配置することを特徴とする電界効果型SiMIS半導体装置の製造方法である。
更に、請求項14によって提供されるのは、Si半導体単結晶基板上に、電界効果型MIS素子を形成するにあたり、該MIS素子のゲート絶縁酸化膜を単結晶とし、しかも、該半導体基板と該単結晶強誘電体膜を含む系に於いて、その系の対象特性を表現するabinitio分子動力学理論に基ずいた評価関数として系全体の自由エネルギを取り上げ、この自由エネルギが、主動作条件下で最小になるように前記単結晶強誘電体膜を配置するとともに、その単結晶強誘電体膜の酸素欠陥濃度を0.01%以下に抑えたことを特徴とする電界効果型SiMIS半導体装置の製造方法である。
更に、請求項15によって提供されるのは、Si半導体単結晶基板上に、電界効果型MIS素子を形成するにあたり、該MIS素子のゲート絶縁酸化膜を単結晶とし、しかも、該半導体基板と該単結晶強誘電体膜を含む系に於いて、その系の対象特性を表現するabinitio分子動力学理論に基ずいた評価関数として系全体の自由エネルギを取り上げ、この系の自由エネルギが、主動作条件下で最小になるように前記単結晶強誘電体膜を配置するとともに、特に局所的なエネルギの凸凹が偏差値3σ以内に納め、その単結晶強誘電体膜の酸素欠陥濃度を0.01%以下に抑えたことを特徴とする電界効果型SiMIS半導体装置の製造方法である。
以下にこれらの最適化手順ならびに、試作結果を記す。
ここでは、単結晶強誘電体膜として、例えばPZT系ペロブスカイトを例に取り、またSi表面を例えば(100)面とした。この基板面は例えば金属電極としてFCC構造のCuでも良い。一番重要な問題は、単結晶強誘電体膜PZT系ペロブスカイトとSi(100)面とをどのような結晶位置関係に配置すれば良いかという問題に置き換えられる。ここでは、本発明の骨子である評価関数として系の全エネルギを取り上げた。しかも、この系としては、Si/PZT界面を取り上げた。本発明者等は厳密な評価関数としてのエネルギ表現式を予め作成しておき、これを用いることにより、印加条件下でのPZTとしてのペロブスカイトと、Si(100)の位置関係は以下の様にするのが最適であることを見い出した。即ち、ペロブスカイトのC4v構造表現のa1軸とa2軸のなす角の中線 を、Si(100)面上の[110]軸方向に合わせ、且つc軸を傾け、第1Siと第3Si位置をSi(100)面上のSiに一致させることが最適値であることを見い出した。またこの時、PZT中には0.01%以内の酸素欠損が許容範囲であることも同時に示した。この様に作成した本発明の提案する概念図を図34に示した。
図34に付いて詳細に説明する。図34aの上半分は、Si側の単一格子の平 面図を示しており、下半分はPZT面のスケッチである。また、図34cは、[ 110]Si方向に対して、ペロブスカイトを結ぶ線を並行に配置させるのが最も良いと指摘している。この時、後で詳細は説明するが、シリコンと結ぶ直線が[110]Siに平行になる為には、ペロブスカイトのC軸が傾く必要がある。
こうすることによって、本発明者らの計算によれば、単に,ペロブスカイトとSiを結ぶ線を並行に配置させて、C軸を傾けない場合にくらべ、全エネルギがさらに10%低下することをも見い出した。この上で更に、実際素子への利用を考え、PZT中に0.01%以内までで有れば、充分単結晶ペロブスカイトの本来の特徴を示す事を見い出した。
ここで指摘する界面エネルギの問題は以下に説明するように、素子特性には多々影響を及ぼすものである。強誘電体膜と基板とが形成する界面に発生する不整合によるエネルギは、不対原子を形成しやすくし、これにより、界面準位の形成を助長することになる。この界面準位は、伝導現象の妨げや、信頼性の低下につながる。また膜中においても、異常ホ゛ント゛等の欠陥ができ、準位になる可能性を はらんでいる。ここでは、これらの欠陥の許容範囲を初めて示すものである。
本発明によって作成した、MIS特性と従来例との比較をも行った。ここでは、特に界面準位を効果的に評価できるものとして、CV法を用いた。本発明による特性には、データのばらつきも正確に示した。本発明によった場合、その界面準位に値は従来例に比べ、約1/10に低下していることがわかる。また本発明 によるMIS素子の移動度も同一電界状況下で評価すると約12%の向上が見られた。また信頼性も向上した。
本発明者等は今までにない新しい軸関係を発明するに至った手順をさらに詳しく説明する。図35が本発明者等が作成した、全く新しい厳密abinitio/分子動力学シミュレータの構造図である。図36の上半分は、abiniti oの部分で、下半分は分子動力学シミュレータの部分である。このシミュレータを用いて、初めて、発明者らは、例えば ペロブスカイトの構造が如何なるものかを示した。図37が、本発明者らが、初めて厳密な構造を示すことができた、ペロブスカイトの透視図である。図中の大きい玉は酸素を示し、小さい玉は、Pbを示している。ここでZ軸はペロブスカイトのC軸を示しており、図でわかる様に、ペロブスカイトのどの面とSiのどの面が巧く接合するかは、この様に複雑な面であるので、極めて推測は困難であることが良く分かる。
また、上記の分子動力学シミュレータの精度チェックとして、酸化膜に着目し、分光データの対比作業を進めた。即ち、時々刻々のSi原子や酸素原子の位置等から、O-Si-OやSi-O-Siの、ロッキンモードやストレッチングモード等を読み とり、これらの数値から固有周波数を求めた。これと実測値とを比べたところ、それぞれ、440cm−1と1100cm−1に対し、490cm−1及び1111cm−1であり、妥当な範囲である。また圧力を印加した場合、O-Si-Oモードの強度が増大することも実測値と良く一致していることが分かった。
これらの検証作業をもすすめながら、本発明者らは、PZTに非常に大きな圧縮応力を印加した場合についても計算してみた。界面としては、この様な構造をとったほうが、全エネルギの利得があることもわかった。
計算に使った領域等について若干詳しく論じてみる。計算ではSi基板側に、深さ方向に数十原子層をとり、また、表面での奥行き方向に数十原子層をとり、また、左右方向にも数十原子層をとり、計算領域をとった。またペロブスカイトにおいても、縦横高さをそれぞれ数十原子層ずつとった。
今までにない手法で、新しい、位置関係を見いだせたのは、この方法で摘めている。本発明者らは、強誘電体膜及びSi単結晶を含む領域に対して、厳密なイオン結合ポテンシャルの厳密積算式を開発した。従来の手法では、厳密性に欠け、特に系のエネルギを正確に求めることが出来なかった。本発明者はこの明細書のなかで従来の欠点や不正確な場所に付いて順次明らかにして行く。
本発明者によると、以下のように説明できる。本発明者らは、従来にない新しいシミュレータを構築し、一例として、強誘電体膜の単結晶、具体的にはペロブスカイトについて計算を行った。
以下に、本発明者が作成した、シミュレータについても従来との比較をしながら、若干説明を加える。たとえば、まず、酸化膜について述べてみる。Si−O間、O−O間、Si−Si間の3種類のポテンシャルを厳密に表現しなければならない。Si−O間とO−O間、Si−Si間のtotalのポテンシャルは実際は すこぶる難しく、3種とも距離rの項により積算量は、無限大に発散することになる。また計算自体もクーロンポテンシャルは遠くまで裾を引くので打ち切ることができず、少し厳密な計算では従来の手法ではEwaldの方法を使っていた。し かし、後述するが、計算に厳密性を欠いている。
実際に配置するに当たり本発明者らが用いた手法を以下に示す。ここでは、本発明者らが初めて作成したシミュレータを以下に示す。
本発明者らは、鋭意検討し、計算物理学に則りながらも、新しい独自のabinitio分子動力学システムを構築した。この新しいシステムがあったからこそ、新しい試みが可能になった訳である。以下に本発明者等が作成したものを詳解する。
ここではSiO2を例に採っているが、他のものについては、単に元素の番号 等を入れ換えるだけで使用できることも確認した。
本発明者等が提案した新しいシミュレータの全体構成を図38に記す。分子動力学部分(MD)と密度汎関数(DFT)とを合成した新しい手法を提案した。
まず本発明者らによる分子動力学法と密度汎関数法との使い分けを図39を用いながら以下に示す。
図36に示す様に、分子動力学(MD)部分はフローチャートの下半分を占め、密度汎関数部分(DFT)は上半分を占めている。分子動力学法は零Kでない所謂有限温度でのやや大きな系を扱い、原子(イオン)間のポテンシャルには、予め第一原理から求めたポテンシャルを用いる。他方、密度汎関数(DFT)は、小さな系の基底状態を求める時に使うことにした。そして、当然DFTの部分は温度は零Kとした。ポテンシャル部分は、電子系の基底状態から求めた。具体的には波動関数(KS方程式)を解き、局所密度汎関数を併用した。
本発明者等によるシミュレータの計算手続きを実際にSiとOを例にとり説明してみる。DFT部分で、縮重のない電子系の基底状態の全エネルギを求めて置
〔外1〕
]である。
ン相互作用エネルギー、第3項は、電子の1体近似を想定して、電子間相互作用のない時の運動エネルギーであり、
と書かれる。第4項は交換相関エネルギーで、その形は局所密度近似による。電
〔外3〕
を自己無撞着(self-consistant) に解かなければならない。
とした。ε
xc(n)の形はWigner他の式など、いろいろ提案されている。さらに、
〔外7〕
のクーロン場をそのまま使うのでなく、核の位置での特異性を消したノルム保存擬ポテンシャルでおきかえる。このようにならすことによって、フーリエ展開での成分の個数をおさえることができた。
全エネルギーとして、本発明者等は厳密を期するため、原子(イオン)間のクーロン・ポテンシャルも含むものを使う。即ち、
ε など)。式(2071)から右辺第1項はDFTの式(2065)から右
〔外9〕
である。その部分はDFTの場合と同様に1体近似の式(2065)からVeff
〔外10〕
ここで新しい考え方を入れた。Eをポテンシャル・エネルギーと見なして、他に仮想的な運動エネルギー
を導入する。従ってラグランジアンは
L=K−E …(2073)
式(2073)中のM
Iは本当の原子(イオン)の質量だがμ、μ は仮想的な 質量であり、後述のように目的に応じて変えた。束縛条件として正規直交条件
を課するから、オイラー方程式を作る時に未定乗数Λ
ikを導入して、Lに
式(2073)の運動エネルギーの温度Tが対応する。古典統計力学のエネルギー等分配則によって各自由度が温度を持っているのだから、Σ1/2MIRI 2との関係でいえば仮想的でなく現実の温度である。またμ、μ の大小によってνiν を抑制したり増強したりできる。
例えば、μ M
Iとして高温からT→0とすれば、R
Iをとめたままψ
iが変化 して電子系の基底状態を得る。この時式(2076)の左辺が0になり、右辺は式(2071)とその下の注意により
となるからψiの線形結合を適当に作る{Λik}が対角線化されてKS方程式( 2064)を得る。つまりDFTでKS方程式を解いたのと同じ結果をsimulated annealing で得たことになる。ただし、温度Tの状態を作るのに、モンテカルロ法によらず、仮想的な力学系の運動を追っているから、dynamical simulated annealing(DSA)と呼ばれる。
式(2076)の積分をする時、必ずしもψiをそのまま変数として扱うので なく、それを平面波で展開した展開係数(つまりフーリエ係数)を変数として扱うことができる。この展開係数の個数をむやみにふやさないためには、DFTの場合と同様に擬ポテンシャルで核の位置での特異性を消しておくことが必要である。特に結晶の周期性を利用できるときには、上記の{Λik}を対角線化するψiの線形結合をφnkと記す。kはブリユアン帯域内の波数ベクトル、nは帯域イ ンデクス。
を満たすはずである。λ
nkは行列{Λ
ik}の固有値であり、エネルギー準位である。φ
nk(n、t)を展開して
式(2075)と式(2076)を式(2074)に代入して
を得る。この式は振動部分の積分を解析的にすましているから、純粋に数値的な方法よりもΔtを大きくして、計算量を減少することが出来る。
クーロンポテンシャルは、以下の様になる。本発明者は分解してゆくと4項に分かれる事を見い出した。特に従来は3項しかなかったが、新しく第4項があることを確認した。これらの方法を順次しるす。これらは以下に示す様に、本発明者らは、まず厳密な式の出発として、誘電率を加味した基本式から解き始めた。
導体ε=∞と真空ε=1の場合とでクーロンポテンシャルが異なり、Lは(立方体の)単位結晶の一稜、Σは単位結晶内でとり、Zi、酸素i第i粒子の荷電と位置である。これは球内の荷電によって球の内面に分極が生じることによる。
導体でない球の内面に双極子の層が出来るが、上記式の最終項がちょうどその効果を打ち消す働きをする。Ewaldの方法は左辺、ε=∞のものを与えるから、真 空内の値を求めるには上記式の最終項を加えなければいけない。結果だけ掲げる。この部分でも従来の導出はしばしば誤りがある。
と表される。ここに、nx,ny,nzがそれぞれx,y,z方向の稜のベクトル で、nx,ny,nzはそれぞれ(バルク結晶の場合)−∞から+∞にわたる整数 である。Σ′の´はn=oの時のj=iを除くことを意味するものとする。
を見い出した。
指数m
x,m
y,m
zはそれぞれ−∞から∞にわたる整数である。 = の時は、
〔外15〕
を使った。
これは を含まないから力には関係ない。さらに次ぎのように変形できる。
今までの結果をまとめて、クーロン・ポテンシャルは、結局、
/Lx,0,0)+my(0,1/Ly,0)+mz(0,0,1/Lz)である。上記の表現の具合のいい点は、もとのΣの項が逆数のオーダーでしか減衰しないのに対して、Φ(1)ではΣの項がerfcの因子によって、Φ(2)ではΣの項がexpの因子によって、急速に減 衰することである。Φ(3)での減衰との遅速に逆向きに効くから、両者のバランスがとれるような適当なκを選ぶ。距離の近いところから順にクーロン力の寄与を計算した結果であり、しかも周囲が導体である場合の結果である。周囲が真空である場合によればもう1項が加わる。これが特に従来省みられなかった部分である。
本発明者等が作成した回転点群を図40を用いて以下に示す。Si原子を半径rの球の中心(0,0,0)に置き、球面上正四面体の頂点の位置に4個の酸素原子を置 く。C原子の位置は、オイラー角(θ、φ)で指定できる。それを(x,y,z)座標で表すには、北極に向かうベクトル(0、0、γ)に対して、まずy軸まわりθの回転、つ いでz軸まわりφの回転をしてやればよい。
ここで、4個の酸素原子のうち1つを北極(0、*)に置いた時、しかもφは不定で*としておくと、他の3つの酸素原子の位置は、未定の角ψによって、それぞれ、(109°28´,ψ)、(109°28´,ψ+120°),(109°28´,ψ+240°)と表される。従って、それらの(x,y,z)座標は、式(2061)で、cosθ=-1/3、sinθ=2√(2)/3、そしてφをそれぞれψ、ψ+120°、ψ+240°としたものを、t(0,0,r)に作用させて、
となる。北極を指していた第1の酸素原子を(θ、φ)の向きにするには、式(2 061)のR(φ、θ)をそのまま作用させればよい。その結果、本発明者等は以 下のように求めた。4個の酸素原子の(x,y,z)座標は、
となる。単位格子内の4個のSi原子の(x,y,z)座標をパラメタa,b,cで表す。図 40と照合すると、2aが縦方向、横方向に共通の周期である。
I ( 0, 0, 0)
II ( a, b, c)
III(a-b,a+b, 2c)
IV ( -b, a, 3c)
V ( 0, 0, 4c)
すなわち、第1原子を原点におき、z座標が一定値cずつ上り、(x,y)面に正 射影すると、I 、II、III、IVで正方形をなしている。第5原子は次の単位格子の第1原子でz座標が4cふえて、第1原子の真上にある。各Si原子に属する4個の 酸素原子の配置は、Si原子に相対的な配置がz軸まわりに90°ずつ回転して行く。それに上述のようなSi原子の移動が加算される。第1Si原子周囲の4個 の酸素原子の(x,y,z)座標は式(2061)をそのまま使えばよく、第2、第3 、第4Si原子の周囲の酸素原子の(x,y,z)座標は、第1原子周囲の酸素原子に 対して、それぞれ
R(0, 90°)の回転 と( a, b, c)の平行移動
R(0,180°)の回転 と(a-b,a+b, 2c)の平行移動
R(0,270°)の回転 と( -b, a, 3c)の平行移動を作用させれば得られる。
その他に、I II1とIII2、II III1とIV2、III IV1とV2の組がそれぞれ同一の酸素原子であるが、これらは(当然のことだが)上と同じ関係式を与える。酸素原子I1とII4は、横方向に1周期分2aだけずれている。
この様に調べてきた上で、さらに実際プロセスを考えてどの程度の欠陥密度まで許容されるかについて調べた。即ち、今、例えば3200個のSiと3200個のTiを用意し、これに0個から100個の酸素欠陥を作った。この場合の時々刻々のSi原子や酸素原子の位置等から、次の事が明らかに成って来た。まず、酸素原子及びTi原子の動的挙動とPZT構造因子との関連を調べてみた。単結晶PZT中に酸素欠陥が存在すると、数原子先の遠方のTi及び酸素までも、無欠陥の場合と比較すると明らかに擾乱を受けており、しかも、点欠陥の移動度はすこぶる大きいことである。これらの動的挙動を、単結晶Siの場合と比べると、PZTの場合の方が極めて影響も範囲も大であることが分かった。この様な活発な動きは、PZTにはイオン結合成分がかなり有り、分極しやすい性格を持っており、従って、局所的な欠損が引き金に成ってポテンシャル変化に連動するものとも考えられる。本発明者らは、このような手続きで、詳細に調べたところ、0.01%以下であれば充分単結晶強誘電体膜の利点を引き出すことが出来ることを見い出した。
本発明者による強誘電体材料の原子レベル計算について以下にします。本発明者等は、PTOやPZO、PZTの原子レベルでの構造把握や、構造と膜特性との把握計算を進めた。計算の結果を報告してみる。原子レベル計算で用いている対象構造としては、(1)PTOの完全結晶で、酸素、Tiがそれぞれ高対称位置 にあるもの、(2)PTOの完全結晶で、酸素が、0.3Aずれたもの、(3)前記(2)に 酸素欠損を、c軸にそった部分に1個入れたもの(図37)、(4)前記(2)に酸素欠損を、a軸にそった部分に1個入れたもの(図43)、(5)PTOのTをZr に置き換えたもの等である。特に、(5)においては、Zrを、人工的に層状に置 換したものも実行した。酸素やPbの働きは、PTOやPZO等で異なっていることが初めて分かった。これらから、強誘電体膜の「自発分極」のきっかけが掴めた。また、強誘電体膜の原子配置や歪み、欠損構造等と、基本膜特性との関連もつかめた。
最も基本的なPbTiO3膜の原子レベルでの計算についてまず示す。引っ張 り歪みを受けた場合や、試みにPb位置を平衡点からずらせた場合(図44、図45)や、酸素欠損を有した場合等について、膜特性を求める上で最も基本になる電荷分布が上手く求まるかを調べた。図44はPbTiO3の単位胞の絵で、 さらにPbを中心位置から、少しずらせた場合である。この時の電荷分布の等濃度図を求めたものである。また図45はやはり単位胞であるが、Ti-酸素結合の 部分をひずませた場合である。この結果から以下の点が分かってきた。(i)全体 に非常にイオン結合成分が多いこと。また、(ii)電気陰性度の強い酸素に殆ど大部分の電荷が引き寄せられ、次にTiの廻りに電荷があることが分かった。また(iii)このままでは電荷分布が対称的で自発分極は起こり得ないが、酸素「欠損」 や、Tiの位置変動を考えると、局所的電気的中性が破られることが、容易にわかる。さらに引っ張りにより、誘電率が低下することも確認出来た。
さらに例えば、25℃で熱擾乱がないとして見たときの計算結果では、以下のことが分かった。即ち、点群C4vからなる理想的なPbTiO3膜の安定構造を求めた。その結果、酸素を極わずかずらしても、エネルギの安定な元の点群C4vに戻ることが分かった。さらに大きく、酸素位置を0.3Åずらしたところにも、 通常言われている様に、もっと安定な平衡点があるか。これを調べた。出発点0.3Åより極わずか離れた位置に向かってさらに安定になりつつある。この辺りが 通常言われている位置で、このずれ故に自発分極が生じるとも言われている。安定点が少なくとも二つあり、1つは単純立方晶の位置にあるもの、もう1つは若干ずれた位置であることが初めて分かった。
分極の起源は、局所的に電荷量の「中性」を、くずす様な原子位置が、自由エネルギ的に安定点であることとも置き換えられる。また、以上の計算では、酸素の存在は極めて大きいことが分かった。それ故、酸素欠損は、ひょっとしたら、結晶的にも、電子論的にも大きなトリガーになりうることがわかる。即ち、
(1)各3種類の原子にそれぞれ電荷分布は局在している様である。
(2)周期律表にある順序に従って、O及びTi及びPbは、それぞれ(1s2,2s2,2p4),(1s2,2s2,2p6,3s2,3p6,3d2),(1s2,2s2,2p6,3s2,3p6,3d10,4s2,4p6,4d10,4f14,5s2,5p6,5d10,6s2,6p2)である。これを反映して、特にPbの電荷分布領域が大きいのが分かる。
(3)良く見ると、とりわけ、Tiの電荷分布は直ぐ隣のOによって変調されて いるのが見える。
これを抑止するため、F,Cl、Br、I等のハロゲン元素を上手く置き換えてやる手も有ることが分かった。
例2 上の例では、Si・PZTについて記したが、単結晶強誘電体膜として、スピネルMgAl2O4膜、酸化セリウムCeO2 、チタン酸ストロンチウムSrTiO3 を用いることもできる。また、酸化アルミニウムAl2O3、酸化ジルコニウムZrO2 、酸化イットリウムY2O3、イットリウム安定化ジルコニウムYSZ、PrO2 、弗化カルシウムCaF2 など、および、その積層膜をも用 いることができる。また、該単結晶強誘電体膜と下地シリコンウエハあるいは下地電極との界面にシリコンPZTを挟む変形は本例以外の構造においても応用することは可能である。
また、電圧下での自由エネルギを小さくすることが大切であり、歪のエネルギだけに着目すれば、そのときは必ずしも最低になっていなくても良い。我々がペロブスカイトについて示したのも、思想に基づいている。
本発明による実際の測定系についても説明する。図47は、ソーヤータウワの回路である。この方法の回路について説明する。Tは、高圧トランス、Coは既 知固定容量、CROは陰極線オシロスコープである。Coを、強誘電体試料の容 量より、充分大きく選ぶと、Tの2次電圧は殆ど試料にかかるが、それは、抵抗Rdによって、適当に分割されて、CROの横軸に加えられるから、CROの横軸の振れは、試料内の電場Eに比例する。一方、CROの縦軸の振れは、Coの 両端子間の電圧を表すが、これは、DS/Coに等しい。即ち、CROの図形は 、D−E曲線を表し、縦軸の振れからD,Eの絶対値も容易に得られる。自発変位Dsは、飽和部分をE=0へ外捜したときのDとして得られ、自発分極Psは、Ps=Dsできまる。損失のない常誘電体をこの回路で観測すると、原点を通る直線になるが、損失のある場合や、強誘電体はヒステレシスをもつ。図48に本発明になるD−E曲線を実線で表し、比較のために従来のPZTの結果を破線で示した。これから分かる様に、本発明になるものは、非常に誘電率が大きいこともわかる。
ここで指摘するところは、強誘電体膜と基板とが形成する界面に発生する不整合によるエネルギを印加時で最低にしようとするものである。しかも、この条件で実際利用上のことを考えその品質の許容範囲も示した。この界面エネルギを縮小することは、伝導現象の妨げになる準位の形成抑止や、信頼性の向上につながる。この構造を本発明によれば、新しいabinitio/分子動力学シミュレ ータによって予測するものである。所望の評価関数を搭載したabinitio/分子動力学シミュレータは、経験的なモデル等を一切用いていないので、この システムに指摘するところは、全て自然現象を完全に予測できるものである。ここでは、系の印加時の全エネルギに着目した例を示しているが、例えば誘電率を所望の値にすべく設計についても、さらには、劣化抑止を所望の問題であれば、それに適した設計が可能になる。
以上2次元の物性を計算する方法を示した。3次元の物性の計算方法の一例については、図8と共に既に説明した。
叉、実際の対比や診断がどの様にして行われているかをしめそう。即ち、シミュレーション結果は、スケジューラ(13)に転送される。この時、他方、実際の実験系からは、モニターデータが、時々刻々送られ来る。そして、これらは、ファイル(4)に格納されるととものに、ファイル(4、3)からのデータ転送は、スケジューラによって同期させておく。そして、対比診断検定部IVにある対比・検定・診断ブロックに転送される。
図49は、本発明に使用されるクラスタ化ツールにおいて、in−situモニタ測定装置がどのようにして付加されているかを示す概念図である。図49に示すように、モニタ量は、例えば、入力ガス構成分析やXPSシグナル、FT−IRシグナル等であり、これらは、所定の時間空間内における平均値であり、ほんの一部分の情報である。本発明者らは、これら限られた測定量と、他方計算から得られた十分な情報とを対比させる方法を提供するものである。これらは、図50の様に、モニター上のグラフィックな画像として表示される。
図51は、本発明によるabinitio分子動力学シミュレーションシステムへの入力データの一例としての温度シーケンスを示すものである。また図52は、本発明によるabinitio分子動力学シミュレーションシステムによる、原子レベル実行予測の概念を示すものである。このように、例えば原子レベルで、基板上での反応や、分子の解離などが記されている。
図53は、本発明によるabinitio分子動力学シミュレーションシステムにより得られた原子レベル実行予測の可視化図である。これは、非晶質Siが堆積されたあと、600℃近傍の高温でアニールされた状態を原子レベルでみたものである。
また、本発明のabinitio分子動力学シミュレーションシステムは、Siのみならず、酸化膜に対しても適用することができる。即ち、図54は、本発明によるabinitio分子動力学シミュレーションシステムのSiO2につ いての実行例を示す。ここでは、SiO2におけるSi−O−Si角が時々刻々 求められており、これから、フーリエ変換し、固有振動数をもとめ、これを重畳することにより、FT−IRのスペクトルが計算のみで求められる。この予測結果を用いて、実測値の時々刻々データと比較検定すれば良い。
図55は、さきに示した図51の入力データを用いた場合の、本発明のシミュレーションシステムを用いた原子レベルでの実行例を示す。図55から、予想どおり実行されていることがわかる。
図56から図69は、図51の入力データを用いた場合のSi膜の成膜から、アニール、さらに固相成長までにおけるSi−Si間距離、結合角を、時間と温度を変化させてもとめたものである。
図1に示したシステムでは、基本膜特性の算出を膜特性算出ユニットIIIで行った。しかし、図70に示したシステムのように、特別な膜特性算出ユニットIIIを省略し、プロセス変動算出ユニットI内部で必要な特性の算出を行ってもよい。それらは設計の要請に応じて、適宜変形可能である。
以上のように、本発明によって、ハードウエアや計算量を従来に比較して大幅に軽減することが可能となった。例えば、計算量について言えば、プロセスシミュレータTOPAZでは、図98に示すように、メモリは約24Mバイトで、実行時間は数分間となる。これに対して、有名なプロセスシミュレータSUPREM4では、粘弾性等を考慮しているので、メモリは200Mバイト以上、実行時間はクレイで1日にも達する。さらに分子動力学やab−initio計算になると、メモリが数倍、実行時間も数倍にも達している。
このような状況では、分子動力学やab−initioツールをフルに用いて、原子の動的挙動の予測や、光学的諸特性の予測には、計算機のコストが大幅にかかかってしまう。本発明者らは、これら不利益を、新しい数学的な技巧を駆使し、克服する方法を見出すとともに、この技巧により、シミュレータの新たな利用技術をも見出したものである。
従来でも、上記に若干の手続きを加えたものもある。それが図99に示すフローシートであり、従来の分子動力学シミュレータにおける密度汎関数の流れの主要部分を示す。即ち、この手法では、仮想質量を与え、全系のラグランジアンを用いる。これは、本発明者等が先に記した式(1223)に該当する。これはすでに知られているCar−Parrinello法に相当するものである。この方法は一見高速に計算できるとされている。しかし、それでも例えば、クレイ社レベルのスーパーコンピュータでも、1つの場合の計算でも1昼夜を要する。従って、種々のパラメータの変動を調べるため、1回1回計算をしていたのでは、到底、実時間での診断ツールとはなり得ない。