JPH11176906A - 電子部品の製造方法、製造システム、設計方法、及び記録媒体 - Google Patents

電子部品の製造方法、製造システム、設計方法、及び記録媒体

Info

Publication number
JPH11176906A
JPH11176906A JP34644097A JP34644097A JPH11176906A JP H11176906 A JPH11176906 A JP H11176906A JP 34644097 A JP34644097 A JP 34644097A JP 34644097 A JP34644097 A JP 34644097A JP H11176906 A JPH11176906 A JP H11176906A
Authority
JP
Japan
Prior art keywords
molecular dynamics
simulation
electronic component
equation
manufacturing
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Withdrawn
Application number
JP34644097A
Other languages
English (en)
Inventor
Shinji Onga
伸二 恩賀
Takako Okada
多佳子 岡田
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Toshiba Corp
Original Assignee
Toshiba Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Toshiba Corp filed Critical Toshiba Corp
Priority to JP34644097A priority Critical patent/JPH11176906A/ja
Publication of JPH11176906A publication Critical patent/JPH11176906A/ja
Withdrawn legal-status Critical Current

Links

Classifications

    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02PCLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
    • Y02P90/00Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
    • Y02P90/30Computing systems specially adapted for manufacturing

Landscapes

  • Container, Conveyance, Adherence, Positioning, Of Wafer (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

(57)【要約】 【課題】 半導体装置製造プロセスを、テストピースな
しに、所望の工程通り又は修正しながら進行することを
可能とする半導体装置の製造方法を提供すること。 【解決手段】 複数の工程からなる半導体装置の製造方
法において、前記複数の工程の少なくとも1つにおける
実観測データを得る工程と、abinitio分子動力
学プロセスシミュレータ又は経験的ポテンシャルを与え
た分子動力学ミュレータにより、前記複数の工程の少な
くとも1つにおける予測データを得る工程と、前記予測
データと実観測データとを逐次、実時間で比較検定する
工程と、前記比較検定により、製造工程因子の設定値
と、前記実観測データから推測される前記複数の製造工
程因子との間に有意差が認められた場合、前記製造工程
因子を逐次実時間で修正処理する工程とを具備すること
を特徴とする。

Description

【発明の詳細な説明】
【0001】
【発明の属する技術分野】更に、本発明は、電子部品の
製造方法、製造システム、設計方法、及び記録媒体に係
り、特に、枚葉式又はクラスタ化した半導体装置製造装
置における限られた量又は質的内容の“その場”測定値
を最大限に拡張し、テストピースを用いることなく、半
導体装置製造プロセスを、所望の工程通りに又は修正し
ながら進行することを可能とする半導体装置の製造方
法、製造システム、設計方法、及び記録媒体に関する。
【0002】
【従来の技術】従来、半導体装置の製造装置は、次のよ
うな変遷を経て現在に至っている。これを各世代の流れ
のなかで見てみると、dRAMの世代とその世代毎のチ
ップコストは、64Kから256Mまで容量が増加する
とともに、図128に示すように変化する。なお、図1
28の縦軸は、1Mにおけるチップコストを1として表
現している。また、図128において、プロットされた
点はチップコストを示し、横軸はdRAMの世代を示し
ている。
【0003】図128に示すように、64Kから256
Kまでは、比較的工程数も安定しており、チップコスト
も安定している。それが、4Mから世代を経る毎に、工
程数が増えており、それに伴ってチップコストも上昇し
ている。この工程数の増大によるコストの増加を少しで
も低減させるため、図128にも示したように、例え
ば、64Mでは、(1)ムダの排除、(2)作業の標準
化、(3)チップサイズの縮小、(4)微細化推進、
(5)生産性向上、(6)FA化推進、(7)大口径化
等の手段を講じて、コスト低減をはかっている。
【0004】256Mになると、新材料の適用や、高誘
電体の採用、新構造の適用、立体化等を図ることによ
り、バイルールに乗せてゆくのである。しかし、1Gク
ラスになると、状況は一転し、バイルールに従うかどう
か不明となる。これまでの256Mクラスまでの製造装
置個々の高性能化を極限まで追求して行く従来の方法で
は、もはや事業としての解を求めることは困難となる。
製造装置個々の高性能化を極限まで追求して行くこと
は、即ちコスト高を招いてしまう。
【0005】また、図129に示すように、各工程毎
に、高性能化された装置が、例えば、工程AからB→C
→D→E→Fの順序で使用されて行くものとする。この
場合、設備領域の増大もコスト高を招くことは明らかで
ある。
【0006】また、1Gクラス以上で、もう一つ忘れて
はならないことは、各工程において物理的必然性による
精密制御が必須になってくることである。例えば、汚染
の遮断、自然酸化膜の排除、原子レベルでの平坦化の追
求等である。現在では、これらのことを考慮して、研究
段階ではあるが新しい概念として、図130に示すよう
なクラスタ化した製造装置が提案されている。このクラ
スタ化した製造装置(クラスタツール)は、バッチ式も
あるが、主として枚葉式であり、また、図130からも
分かるように、例えば中央に搬送系があり、各工程は、
それぞれAからB→C→Dの如く、矢印に示すように順
次進行して行くものである。このとき、各工程及び各工
程内では、多くの場合高真空に保たれ、自然酸化膜の被
着も減少でき、また、汚染の遮断、原子レベルでの表面
制御ができるようになってきた。
【0007】一方、半導体上に形成した電界効果MIS
素子について、昨今の高速化及び高集積化の技術開発の
中で、勢い、強誘電体は薄膜化に向かっている。強誘電
体膜を薄膜化すると、勿論閾値電圧は浅くなり、この分
だけ主動作速度は早くなり、これにより特にAC特性は
格段に改良させる。また、EEPROM等を考えてみる
と、微細化されてくると、非常に過酷な使用状況になっ
てくる。このような場合には従来の通常の製法の強誘電
体膜では、もはや、充分な信頼性がえられない。
【0008】しかし、MIS素子が微細化されても、強
誘電体膜自体の特性の改善策は余り省みられず、電源電
圧がさほど下がらない状態では、特に主動作時にゲート
強誘電体膜に高電界が掛かる。また、チャネル領域から
インパクトイオン化等で発生した電子正孔は、それぞれ
ゲート電極の極性及びドレイン電圧等の境界条件によ
り、強誘電体膜中に注入される。そしてこれらのキャリ
アは、強誘電体膜中にトラップされ長期信頼性を低下さ
せるだけでなく、耐圧低下等を引き起こす。
【0009】もう少し原子レベルでみてみる。熱強誘電
体膜においては、例えばシリコン強誘電体膜に高電界を
印加すると、ネットワークを構成するSi−O結合が、
外部から印加された高電界と相互作用する。そして、結
合が切断されて、電子や正孔が捕獲される捕獲中心が形
成され、続いて通過する電子や正孔は、この捕獲中心に
捕獲されて、膜厚方向の電界強度分布が局部的に平均電
界より高くなり、やがて絶縁破壊に至る。
【0010】これらの状況を改善しようとして、最近で
は、強誘電体膜を単結晶化するアイデアが学会レベルで
はあるが提唱され出した。例えば、シリコン(111)
面上には、酸化セリウムCeO2が単結晶成長すること
が、J.Appl.Phys.,Aol.69(12),p8313(1991)に報告され
ている。あるいは、弗化カルシウムCaF2がシリコン
単結晶上に単結晶成長することが、Japan.J.Appl.Phys.
Suppl.,Aol.21-1,p187(1982)に報告されている。これら
は、一部であり、まだアイデアの状態で有るばかりでは
なく、その指針のもとになる計算手法にも多々疑問が残
っている。それで、一つには、例えば強誘電体膜を例に
上げても、それを用いた強誘電体膜の単結晶の姿自体が
充分認識されていない状況にある。この単結晶強誘電体
膜の構造に付いて報告したものを見てみよう。一つに
は、A. Miyamot, K.Takeichi, T.Hattori, M.Kubo, and
T.Inui, "Mechanism of Layer-Layer Homoepitaxial
Growth of SrTiO3(100) as Investigated by Molecular
Dynamics and Computer Graphics". Jpn. J. Appl. Ph
ys. Vol 31, (1992) 4463-4464.がある。ここでは、ペ
ロブスカイトを単結晶強誘電体膜として、それを下地S
i基板と接着させるとき、安定構造を計算している。し
かし、Si−O−Ti叉は、Ti−Si−Sr間の角度
や距離の初期配置を誤って用いている。また、これら
は、単結晶強誘電体膜の構造を単純化しており、必ずし
もその指針は正しくはない。
【0011】これらのことから分かる様に単結晶強誘電
体膜をどの様に設計するかに至っては殆ど正しい議論で
きていないのが現状である。また他方、強誘電体膜その
ものについては、その特性の向上に関しては、形成プロ
セス上のリファインに未だ終始しているのが現状であ
る。例えば、出来るだけ清浄な面を予め用意するとか、
或いは、膜生成の為のスパッタ温度やスパッタ雰囲気の
議論に終始している。これらは、単結晶強誘電体膜や強
誘電体膜に関する構造を含めての認識がまだ極めて低い
と言わざるを得ない。また強誘電体の自発分極発現の理
由や、その結晶構造との関連が分かっていないので、現
象論的追求にしか過ぎない。
【0012】
【発明が解決しようとする課題】しかし、以上のような
新しいクラスタツールにおいて、現在直面している最も
大きな困難は、その状況観察手法である。即ち、単一行
程を行う装置の場合では、まずテストピースと称するウ
エハを、各工程に導入していた。このテストピースを用
いて、たとえば、膜厚や、屈折率等を測定し、所定の条
件通り工程が進行しているかどうかをテストしていた。
【0013】かかるテストピースの使用による各クラス
タツール内の条件のモニタリングには、量的、質的理解
上において問題がある。例えば、酸化膜の成膜工程又は
成長工程においては、XPSやFT−IRの装置を用い
るのが通例である。これをin−situで測定してい
たとしても、そのデ−タの意味するところは、現状では
十分に検討出来ていない。せいぜいそのスペクトルのピ
ーク波数位置から判断して、所定の物質が生成している
かどうかという判定程度である。このような状況では、
研究者の判断に頼っているところが多い。また、定量的
な対応が出来ず、更に、迅速な対応も出来ない。また、
研究者の経験に頼ることのできない新しい未知の物質が
あるときは、全く判断の仕様がない。
【0014】又、一方で分子動力学シミュレータの利用
も提案されている。しかし、一般に、分子動力学シミュ
レータは、プロセス変動算出プログラムに比べ、(1)
非常に計算が遅く、また(2)対象とする物質の実行領
域が、非常に小さいと言われている。それで、分子動力
学シミュレータと、プロセス変動算出プログラムを直接
つなげると言うことは、殆ど考えられないことであっ
た。即ち、分子動力学部分は、極めて計算が遅いので、
合体化させても、結局は、分子動力学部分の遅いスピー
ドに全体が引っ張れ、実用に向かないと思われていた。
また、分子動力学部分では、計算領域が極めて小さく、
せいぜい数十オングストローム程度であり、一方、プロ
セス変動算出プログラムでは、数ミクロン領域を計算し
ているため、計算領域のずれの点からもこれらを直接つ
なげることは難しい。
【0015】たとえば、本発明者等は、単結晶強誘電体
膜の設置条件として、界面を含む系の全体の自由エネル
ギに着目し、印加条件下でこのエネルギが、最小になる
ように、位置を決定し、これにより、高信頼性のMIS
素子を作るとともに、界面に準位を形成しないようにし
た。特に、系の全体のエネルギを求めるにあたり、単結
晶強誘電体膜の構造を確実に再現することと、これらの
構成している各原子がどの様に運動するのかを厳密にし
かも正確に掴む必要があったが、本発明者等は、ab-ini
tio分子動力学を完成させ、さらに、分子動力学におい
ては独自のポテンシャル積算法を確立し、初めてこれら
の計算を可能にした。しかも、これを用いて、電気的特
性等の向上点を十分保証する意味から、単結晶強誘電体
膜の品質の度合いまでを明らかにしようと言うものであ
る。
【0016】例えば、トンネル絶縁膜に、高電界を印加
すると、トンネル絶縁膜の基本構造を構成している結合
(例えば、シリコン強誘電体膜の場合にはSi−O結
合)と相互作用し、結合が切断されて、切断された結合
に、続いてトンネル絶縁膜を通過する電子や正孔が捕獲
されて、膜中の電荷分布に応じて、局所的に電界が平均
電界より大きくなり、劣化が更に進行する。劣化の進行
を抑えるには、本発明者によれば、この電界による絶縁
膜の結合(例えばSi−O結合)の切断を抑えることで
ある。
【0017】電界と絶縁膜の結合との相互作用は、電界
の向きと絶縁膜の結合の向きに依存する。したがって、
相互作用を弱くするには、絶縁膜の結合の向きが相互作
用が強い向きのものをより少なくなるようにすることで
ある。非晶質のように、全くランダムとした場合、絶縁
膜の結合の一部は必ず電界との相互作用が強い向きとな
る。
【0018】本発明の目的は、強誘電体絶縁膜の結合の
向きを制御し、電界との相互作用を抑え、電界による誘
電性の劣化の抑制を図ることを目的としている。
【0019】これは、高信頼誘電体膜作成技術に関する
もので、本発明者等が、新しいabinitio分子動
力学システムによる計算結果と独自の実験データを基に
生み出したものである。単結晶強誘電体膜は、これを強
誘電体膜等に適用した場合、材料そのものの選択の他、
基板や基板電極との位置関係等を最適化しないと、かえ
って非晶質誘電体膜よりも悪くなり得る。この時の材料
及び結晶方位のクライテリアを示そうと言うものであ
る。本発明者による新しい指針としては、あえて一言で
記せば「主印加電界下で系の自由エネルギが最小になる
様に結晶配置を行うこと等」である。また実際上の問題
として「単結晶強誘電体膜に存在し得る結晶欠陥濃度の
許容範囲」等をも記している。さらに、「主印加時の系
の自由エネルギが最低になる配置のほかに、劣化のトリ
ガを抑止するため、局所的にみた場合の自由エネルギ分
布の凸凹の許容範囲等」をも記している。
【0020】単結晶強誘電体膜一つ取り上げても、4回
対称のペロブスカイトなど非常に複雑な構造を持ってい
ることや、基板との接合面の組み合わせによって大きく
特性がことなる事など、大変複雑であるが、本発明者等
によれば、結晶学的には「損」をしても主動作電圧下で
「得」をしようとするである。
【0021】又、本発明は、クラスタ化におけるモニタ
機能の問題点を解決し、半導体装置製造装置、特に枚葉
式あるいはクラスタ化製造装置における限られた量又は
質的内容の“その場”測定値を最大限に拡張し、半導体
装置製造プロセスを、テストピースなしに、所望の工程
通り又は修正しながら進行することを可能とする半導体
装置の製造方法、製造装置、シミュレーション方法、及
びシミュレータを提供することを目的とする。
【0022】
【課題を解決するための手段】本発明者らは、abin
itio分子動力学を用い、かつ新たな変分計算法を開
発し、(i)各プロセス要素の揺らぎや、各要素のずれ
の重畳現象も診断し、(ii)所望のシーケンス通り実行
されているかのチェックや、(iii)変動因子の発見も可
能とするシミュレーションシステムを開発した。かかる
システムのさらなる特徴は、(vi)プロセスの変動を数
学的・統計的に、しかも効率よく取り扱えるようにした
点と、(v)新しい推論エンジンを付加した点と、(v
i)新しい逆方向システムを付加した点にある。
【0023】更に、本発明(請求項1)は、電子部品の
製造プロセスの製造装置内物理現象のシミュレーション
方法において、前記製造装置内に形成された電子材料の
特性に関する製造工程因子を入力値とし、この入力値を
基に3次元のシミュレーションを行ない製造装置内に形
成された電子材料の特性の予測データを得る方法であっ
て、前記シミュレーションは、分子動力学に基づいて行
われる工程と、流体モデルに基づいて行われる工程とか
らなり、前記分子動力学に基づいて行われるシミュレー
ション結果として出力される3次元物理量と前記流体モ
デルに基づいて行われるシミュレーション結果として出
力される3次元物理量の少なくとも一方は、他方に転送
されそこでのシミュレーションに利用されることを特徴
とする電子部品の製造プロセスのシミュレーション方法
を提供する。
【0024】又、本発明(請求項2)は、上記請求項1
において、前記3次元分子動力学シミュレーションから
前記3次元流体モデルのプロセスシミュレーションへ、
粘弾性係数、応力定数、拡散定数、Cij等の物理量が転
送されることを特徴とする電子部品製造プロセスシミュ
レーション方法を提供する。
【0025】更に、本発明(請求項3)は、上記請求項
1において、前記3次元流体モデルのプロセスシミュレ
ーションから前記3次元分子動力学シミュレーション
へ、粘弾性係数、応力定数、拡散定数、Cij等の物理量
が転送されることを特徴とする電子部品製造プロセスシ
ミュレーション方法を提供する。
【0026】更に、本発明(請求項4)は、上記請求項
1において、前記電子部品は、半導体装置であることを
特徴とする電子部品製造プロセスシミュレーション方法
を提供する。
【0027】更に、本発明(請求項5)は、上記請求項
1において、更に、前記設定値のゆらぎを入力する段階
を備え、前記分子動力学シミュレーションは、前記入力
値とそのゆらぎを基に変分計算を行うことを特徴とする
電子部品製造プロセスシミュレーション方法を提供す
る。
【0028】更に、本発明(請求項6)は、上記請求項
1において、前記分子動力学シミュレーションは、分子
動力学法と密度汎関数法との組み合わせで行われること
を特徴とする電子部品製造プロセスシミュレーション方
法を提供する。
【0029】更に、本発明(請求項7)は、上記請求項
6において、前記分子動力学法は、高階巡回偏微分法に
基づいていることを特徴とする電子部品製造プロセスシ
ミュレーション方法を提供する。
【0030】更に、本発明(請求項8)は、上記請求項
7において、前記高階巡回偏微分法は、3体問題を扱っ
ていることを特徴とする電子部品製造プロセスシミュレ
ーション方法を提供する。
【0031】更に、本発明(請求項9)は、上記請求項
1において、前記分子動力学シミュレーションは、経験
的ポテンシャルを用いて行われることを特徴とする電子
部品製造プロセスシミュレーション方法を提供する。
【0032】更に、本発明(請求項10)は、上記請求
項9において、前記経験的ポテンシャルはATTAポテ
ンシャルであるすることを特徴とする電子部品製造プロ
セスシミュレーション方法を提供する。
【0033】更に、本発明(請求項11)は、上記請求
項9において、前記経験的ポテンシャルはBKSポテン
シャルであるすることを特徴とする電子部品製造プロセ
スシミュレーション方法を提供する。
【0034】即ち、上記請求項1乃至請求項11に記載
の電子部品製造プロセスシミュレーション方法によれ
ば、流体モデルに基づいて行われるプロセスシミュレー
ションによって数ミクロンレベルの広領域の工程誘起応
力の計算を結果を、分子動力学に基づいて行われるプロ
セスシミュレーションによってより細かい部分の検証を
行うことにより、改めて原子の動きを確認する。又、更
に必要に応じて、再び流体モデルのプロセスシミュレー
ションを行い、これを繰り返すことにより、より正確で
より多くの情報を得ることができる。
【0035】又、本発明(請求項12)は、電子部品の
製造工程において、製造装置内に形成された電子材料の
特性を測定して実観測データを得る工程と、前記製造装
置内に形成された電子材料の特性に関する製造工程因子
を入力値とし、この入力値を基に3次元のシミュレーシ
ョンを行ない前記製造装置内に形成された電子材料の特
性の予測データを得る工程と、前記予測データと実観測
データとを比較検定する工程と、前記比較検定により、
前記製造工程因子の設定値と、前記実観測データから推
測される前記製造工程因子との間に有意差が認められた
場合、前記製造工程因子を修正処理する工程とを備え、
前記シミュレーションは、分子動力学に基づいて行われ
る工程と、流体モデルに基づいて行われる工程とからな
り、前記分子動力学に基づいて行われるシミュレーショ
ン結果として出力される3次元物理量と前記流体モデル
に基づいて行われるシミュレーション結果として出力さ
れる3次元物理量の少なくとも一方は、他方に転送され
そこでのシミュレーションに利用されることを特徴とす
る電子部品の電子部品の製造方法を提供する。
【0036】更に、本発明(請求項13)は、上記請求
項12において、前記実観測データから推測される前記
製造工程因子との間の有意差が許容範囲を越えていて、
修正処理が不能である場合には、不良品として処理する
ことを特徴とする電子部品製造プロセスシミュレーショ
ン方法を提供する。
【0037】更に、本発明(請求項14)は、上記請求
項12において、前記製造装置はクラスタ化しているこ
とを特徴とする電子部品の製造方法を提供する。
【0038】更に、本発明(請求項15)は、上記請求
項14において、前記クラスタ化した製造装置は枚葉式
であることを特徴とする電子部品の製造方法を提供す
る。
【0039】更に、本発明(請求項16)は、上記請求
項12において、前記3次元分子動力学シミュレーショ
ンから前記3次元流体モデルのプロセスシミュレーショ
ンへ、粘弾性係数、応力定数、拡散定数、Cij等の物理
量が転送されることを特徴とする電子部品製造プロセス
シミュレーション方法を提供する。
【0040】更に、本発明(請求項17)は、上記請求
項12において、前記3次元流体モデルのプロセスシミ
ュレーションから前記3次元分子動力学シミュレーショ
ンへ、粘弾性係数、応力定数、拡散定数、Cij等の物理
量が転送されることを特徴とする電子部品製造プロセス
シミュレーション方法を提供する。
【0041】更に、本発明(請求項18)は、上記請求
項12において、前記特性は、半導体基板上に形成され
た膜、半導体基板表面、及び半導体基板内部の少なくと
も1つに関する特性であることを特徴とする電子部品の
製造方法を提供する。
【0042】更に、本発明(請求項19)は、上記請求
項12において、前記予測データと実観測データとの比
較検定は逐次実時間で行われることを特徴とする電子部
品の製造方法を提供する。
【0043】更に、本発明(請求項20)は、上記請求
項12において、前記製造工程因子を修正処理する工程
は逐次実時間で行われることを特徴とする電子部品の製
造方法を提供する。
【0044】更に、本発明(請求項21)は、上記請求
項12において、前記分子動力学シミュレーションは、
分子動力学法と密度汎関数法との組み合わせで行われる
ことを特徴とする電子部品の製造方法を提供する。
【0045】更に、本発明(請求項22)は、上記請求
項21において、前記分子動力学法は、高階巡回偏微分
法に基づいていることを特徴とする電子部品の製造方法
を提供する。
【0046】更に、本発明(請求項23)は、上記請求
項22において、前記高階巡回偏微分法は、3体問題を
扱っていることを特徴とする電子部品の製造方法を提供
する。
【0047】更に、本発明(請求項24)は、上記請求
項12において、前記分子動力学シミュレーションは、
経験的ポテンシャルを用いて行われることを特徴とする
電子部品の製造方法を提供する。
【0048】更に、本発明(請求項25)は、上記請求
項24において、前記経験的ポテンシャルはATTAポ
テンシャルであるすることを特徴とする電子部品の製造
方法を提供する。
【0049】更に、本発明(請求項26)は、上記請求
項24において、前記経験的ポテンシャルはBKSポテ
ンシャルであるすることを特徴とする電子部品の製造方
法を提供する。
【0050】更に、本発明(請求項27)は、上記請求
項12において、前記分子動力学シミュレーションで得
られた前記複数の特性の少なくとも1つの予測データを
用いて、更に流体モデルのプロセスシミュレーションが
行われることを特徴とする電子部品の製造方法を提供す
る。
【0051】更に、本発明(請求項28)は、上記請求
項27において、前記分子動力学シミュレーションと流
体モデルのプロセスシミュレーションは、共に3次元シ
ミュレーションであることを特徴とする電子部品の製造
方法を提供する。
【0052】更に、本発明(請求項29)は、上記請求
項27において、前記3次元分子動力学シミュレーショ
ンから前記3次元流体モデルのプロセスシミュレーショ
ンへ、粘弾性係数、応力定数、拡散定数、Cij等の物理
量が転送されることを特徴とする電子部品の製造方法を
提供する。
【0053】即ち、上記請求項12乃至請求項29に記
載の電子部品の製造方法によれば、流体モデルに基づい
て行われるプロセスシミュレーションによって数ミクロ
ンレベルの広領域の工程誘起応力の計算を結果を、分子
動力学に基づいて行われるプロセスシミュレーションに
よってより細かい部分の検証を行うことにより物性定数
の修正量を計算しながら、改めて、原子の動きを確認す
る。従って、実時間での診断を行いながら、所望のシー
ケンス通り実行されているかのチェックや、変動因子の
発見も可能となる。
【0054】又、本発明(請求項30)は、電子部品の
製造装置と、前記製造装置内に形成された前記電子部品
の電子材料の特性を測定する測定装置と、前記電子部品
の製造工程因子の設定値を入力値とし、この入力値を基
にシミュレーションを行ない前記特性の予測データを得
るコンピュータシステムと、前記測定装置及び前記製造
装置と前記コンピュータシステムを接続する手段とを備
え、前記コンピュータシステムは、前記予測データと実
観測データとを比較検定を行ない、前記比較検定によ
り、前記製造工程因子の設定値と、前記実観測データか
ら推測される前記製造工程因子との間に有意差が認めら
れた場合、前記製造装置の制御手段を介して製造工程因
子を修正処理することを特徴とする電子部品の製造シス
テムを提供する。
【0055】更に、本発明(請求項31)は、上記請求
項30において、前記製造装置はクラスタ化しているこ
とを特徴とする電子部品の製造システムを提供する。
【0056】更に、本発明(請求項32)は、上記請求
項31において、前記クラスタ化した製造装置は枚葉式
であることを特徴とする電子部品の製造システムを提供
する。
【0057】更に、本発明(請求項33)は、上記請求
項30において、前記シミュレーションは、分子動力学
法と密度汎関数法との組み合わせによる分子動力学シミ
ュレーションに基づいて行われることを特徴とする電子
部品の製造システムを提供する。
【0058】更に、本発明(請求項34)は、上記請求
項30において、前記特性は、半導体基板上に形成され
た膜、半導体基板表面、及び半導体基板内部の少なくと
も1つに関する特性であることを特徴とする電子部品の
製造システムを提供する。
【0059】更に、本発明(請求項35)は、上記請求
項30において、前記予測データと実観測データとの比
較検定は逐次実時間で行われることを特徴とする電子部
品の製造システムを提供する。
【0060】更に、本発明(請求項36)は、上記請求
項30において、前記製造工程因子を修正処理する工程
は逐次実時間で行われることを特徴とする電子部品の製
造システムを提供する。
【0061】更に、本発明(請求項37)は、上記請求
項30において、更に、前記入力値として前記設定値の
ゆらぎが入力され、前記分子動力学シミュレーション
は、前記入力値とそのゆらぎを基に変分計算を行うこと
を特徴とする電子部品の製造システムを提供する。
【0062】更に、本発明(請求項38)は、上記請求
項30において、前記分子動力学法は、高階巡回偏微分
法に基づいていることを特徴とする電子部品の製造シス
テムを提供する。
【0063】更に、本発明(請求項39)は、上記請求
項38において、前記高階巡回偏微分法は、3体問題を
扱っていることを特徴とする電子部品の製造システムを
提供する。
【0064】更に、本発明(請求項40)は、上記請求
項30において、前記分子動力学シミュレーションは、
経験的ポテンシャルを用いた方法を併用することを特徴
とする電子部品の製造システムを提供する。
【0065】更に、本発明(請求項41)は、上記請求
項40において、前記経験的ポテンシャルはATTAポ
テンシャルであるすることを特徴とする電子部品の製造
システムを提供する。
【0066】更に、本発明(請求項42)は、上記請求
項40において、前記経験的ポテンシャルはBKSポテ
ンシャルであるすることを特徴とする電子部品の製造シ
ステムを提供する。
【0067】更に、本発明(請求項43)は、上記請求
項30において、前記分子動力学シミュレーションで得
られた前記特性の少なくとも1つの予測データを用い
て、更に流体モデルのプロセスシミュレーションが行わ
れることを特徴とする電子部品の製造システムを提供す
る。
【0068】更に、本発明(請求項44)は、上記請求
項43において、前記分子動力学シミュレーションと流
体モデルのプロセスシミュレーションは、共に3次元シ
ミュレーションであることを特徴とする電子部品の製造
システムを提供する。
【0069】更に、本発明(請求項45)は、上記請求
項43において、前記3次元分子動力学シミュレーショ
ンから前記3次元流体モデルのプロセスシミュレーショ
ンへ、粘弾性係数、応力定数、拡散定数、Cij等の物理
量が転送されることを特徴とする電子部品の製造システ
ムを提供する。
【0070】即ち、上記請求項30乃至請求項45に記
載の電子部品の製造システムによれば、コンピュータに
よるシミュレーションによる実時間での診断ツールとし
て、様々な物理問題をリアルタイムで予測することが可
能となる。
【0071】又、本発明(請求項46)は、電子部品の
製造プロセスの製造装置内物理現象のシミュレーション
を行うプログラムであって、前記製造装置内に形成され
た電子材料の特性に関する製造工程因子の設定値を入力
する手順と、この入力値を基に分子動力学法と密度汎関
数法との組み合わせによる分子動力学シミュレーション
を行い、前記特性の予測データを得る手順とを実行させ
るプログラムを記録したコンピュータ読み取り可能な記
録媒体を提供する。
【0072】更に、本発明(請求項47)は、上記請求
項46において、前記分子動力学シミュレーションは、
分子動力学法と密度汎関数法との組み合わせによること
を特徴とする記録媒体を提供する。
【0073】更に、本発明(請求項48)は、上記請求
項47において、前記分子動力学法は、高階巡回偏微分
法に基づいていることを特徴とする記録媒体を提供す
る。
【0074】更に、本発明(請求項49)は、上記請求
項48において、前記高階巡回偏微分法は、3体問題を
扱っていることを特徴とする記録媒体を提供する。
【0075】更に、本発明(請求項50)は、上記請求
項47において、前記電子部品は、半導体装置であるこ
とを特徴とする記録媒体を提供する。
【0076】更に、本発明(請求項51)は、上記請求
項47において、更に、前記設定値のゆらぎを入力する
手順を備え、前記分子動力学シミュレーションは、前記
入力値とそのゆらぎを基に変分計算を行うことを特徴と
する記録媒体を提供する。
【0077】更に、本発明(請求項52)は、上記請求
項47において、前記分子動力学シミュレーションは、
経験的ポテンシャルを用いた方法を併用することを特徴
とする記録媒体を提供する。
【0078】更に、本発明(請求項53)は、上記請求
項52において、前記経験的ポテンシャルはATTAポ
テンシャルであるすることを特徴とする記録媒体を提供
する。
【0079】更に、本発明(請求項54)は、上記請求
項52において、前記経験的ポテンシャルはBKSポテ
ンシャルであるすることを特徴とする記録媒体を提供す
る。
【0080】更に、本発明(請求項55)は、上記請求
項47において、前記分子動力学シミュレーションで得
られた前記特性の予測データを用いて、更に流体モデル
のプロセスシミュレーションが行われることを特徴とす
る記録媒体を提供する。
【0081】更に、本発明(請求項56)は、上記請求
項55において、前記分子動力学シミュレーションと流
体モデルのプロセスシミュレーションは、共に3次元シ
ミュレーションであることを特徴とする記録媒体を提供
する。
【0082】更に、本発明(請求項57)は、上記請求
項55において、前記3次元分子動力学シミュレーショ
ンから前記3次元流体モデルのプロセスシミュレーショ
ンへ、粘弾性係数、応力定数、拡散定数、Cij等の物理
量が転送されることを特徴とする記録媒体を提供する。
【0083】即ち、上記請求項46乃至請求項57に記
載の記録媒体によれば、プロセス装置を実際に稼働させ
ることも、特別な準備や費用を書けることもなく、コン
ピュータ内部に仮想的な半導体製造システムを構築し、
プロセスシミュレーションあるいはを行なうことができ
る。例えば、様々な分子構造モデルを用いて思考実験を
行ない、それを用いたデバイスの特性を予測することに
よって、新たな材料の探索を手軽でしかも効果的に行な
うことを可能となる。
【0084】本発明(請求項9)は、abinitio
分子動力学プロセスシミュレータ又は経験的ポテンシャ
ルを与えた分子動力学シミュレータにより、半導体装置
の製造プロセスの複数の工程の少なくとも1つにおけ
る、複数の製造工程因子の少なくとも1つの設定値、及
びこの設定値のゆらぎを入力値とし、この入力値を基に
変分計算を行い、前記複数の工程の少なくとも1つにお
ける製造装置内に形成された膜、半導体基板上に形成さ
れた膜、半導体基板表面、及び半導体基板内部の少なく
とも1つに関する複数の特性の少なくとも1つの予測デ
ータを得る手段と、この予測データと前記複数の工程の
少なくとも1つにおける前記複数の特性の少なくとも1
つを測定して得られた実観測データとを逐次、実時間で
比較検定して、前記製造工程因子の診断及び判定を行う
手段とを具備する、半導体装置の製造工程因子を診断・
判定するシミュレータを提供する。
【0085】一方、本発明者等は分子動力学シミュレー
タについて、前記したような非常に計算が遅いという問
題、及び対象とする物質の実行領域が、非常に小さいと
いう問題に対してそれぞれ新しい考えを導入して克服
し、分子動力学シミュレータとプロセス変動算出プログ
ラムとを新しい方法で合体させ、しかも、この合体によ
り、今までにない新しい機能を引き出すことが出来た。
以下に順に説明する。
【0086】まず(1)の計算速度の問題に関して少し
詳しく触れてみる。分子動力学では、原子一個一個につ
いて、いわばニュートン力学に従って、運動方程式を立
てて、これを、着目している100個または1000個
の原子全部に渡って連立させて解くものである。
【0087】この運動方程式の基になる力の計算が一番
複雑である。これが、従来の計算手法では遅かったゆえ
んである。各原子は、その廻りの原子から力を受けてい
る。近くの原子からも、遠くの原子からも、力を貰って
いる。これを考慮する必要がある。一般にはテレソフの
式に従って導出されるポテンシャルエネルギVijをもと
め、これを基に、力を算出して行く。しかし、従来で
は、概ね、ポテンシャルを、距離や方向等をもとに、予
め各原子に加わる力について方向毎のテーブルを作って
おく手法を取っていた。この方法では、上記の100個
または1000個の原子について、毎回、テーブルを参
照し、ポテンシャルの大きさを見て行くものである。こ
れは、非常に時間の掛かるものである。またポテンシャ
ルそのものからは、力を算出することは出来なくて、ポ
テンシャルの大きさのそれぞれの方向に対する傾きが、
それぞれの方向の力に成るわけである(図2参照)。こ
れゆえ、非常に計算機の時間がかかることが分かる。
【0088】なぜこの様に、数値テーブルを作るかと言
うと、ポテンシャルの形状を表した式が非常に複雑であ
るからである。上にも述べた様に、ポテンシャル(V)
は、廻りの原子の位置はお互いの角度(θ)、さらに
は、距離(r)等が複雑にからみ合っている。これ故
に、これを、数学的に、微分して関数を作るのは非常に
難しい。例えばθはrの関数であるのに、また廻り回っ
て再びrの関数に成ったりする。
【0089】本発明者らは、この様な一般に非常に複雑
で、上記の様にポテンシャルの表現式の引数が、巡回す
る様な場合について、従来、偏微分が試みられなかった
部分において、新しい高階巡回偏微分式を作成した。こ
れにより、式自体の変形は一見非常に複雑になるが、計
算は非常にはやく出来るようになった。従来のように、
数表を参照しなくて良いようになった。
【0090】では、この高階巡回偏微分式の一番の概念
を以下に示す。これは、一般に3個の原子を想定して、
考えている。すなわち、原子3個に関して一般的に記述
すれば、あとは、この手法を広げて行けばよい。3個の
原子があれば、各原子との距離や、その着目している原
子の間の角度等も算出できる。これらにより、ポテンシ
ャルを一義的に求められ、これをもとに、何れの方向に
も、ポテンシャルの「傾き」を容易に算出することが出
来る。本発明者等は、巡回して用いている変数について
全て偏微分式を作っておき、これらを、連立させて解く
ことで、正しく、微分型が作られていることを確認し
た。
【0091】以下本発明によるシステムの概要を具体的
に説明する。最初に、図1を参照して、本発明の応用の
実際を説明する。即ち、プロセス装置Vは、実際に稼働
している半導体製造システムの一部である。ここでは、
スペクトルメータVIによって光の吸収スペクトルをみ
ている。
【0092】スペクトルメータVIによって得られた光
の吸収スペクトルは、一旦ファイル4として格納され、
スケジューラVIIに送られる。スケジューラVIIに
は、後述する本発明による予測値も入力され、これらは
対比診断検定部IVで比較される。
【0093】対比診断検定部IVでの比較結果が、許容
値を超えていれば、スペックを満たさないものとしてラ
インから外される。また、満たしていれば、その変動分
を推論エンジンとしてのプロセス変動算出ユニットIに
送られる。
【0094】プロセス変動算出ユニット1では、基本的
には原子レベルまで考慮しないマクロ的な性質を計算
し、ファイル2に格納する。ファイル2の内容は、分子
動力学シミュレータIIで原子間位置ベクトルと結合角
が算出される。
【0095】分子動力学シミュレータIIで算出された
原子間位置ベクトルと結合角は、膜特性算出ユニットI
IIに入力され、ここで基本膜特性が算出されてファイ
ル3に戻される。膜特性算出ユニットIIIで算出され
た基本膜特性は、プロセス変動算出ユニットIにも戻さ
れ、パラメータの一部とされる。
【0096】このように、本発明によって、半導体装置
製造プロセスを、テストピースなしに、所望の工程どお
り、または修正しながら進行する事を可能とする半導体
装置の製造方法が実現される。
【0097】実際には、分子動力学シミュレータIIの
内容が特に重要である。具体的には、図7のブロック図
に示したように、2つの部分からなっている。第1の部
分IIaは、本件発明者によって導入された高階巡回偏
微分法を利用するものであり、第2の部分IIbは、三
次以上の高次項を考慮した解析を行うものである。な
お、以後、高階巡回偏微分法や三次以上の高次項を考慮
した解析の詳細が順次説明される。
【0098】即ち、本発明者らは、純粋数学理論にたち
帰り、新たな変分計算法を確立するとともに、新しい計
算化学理論により、各種スペクトルをも予測する技術を
開発した。以下、本発明における新しいabiniti
o分子運動学について説明する。
【0099】本発明者等が開発した新しいシミュレーシ
ョンの全体構成を図3に示す。分子動力学部分(MD)
と密度汎関数(DFT)とを合成した新しい手法を提案
した。まず本発明者らによる分子動力学法と密度汎関数
法との使い分けを図3を用いながら以下に示す。
【0100】図3に示すように、分子動力学(MD)部
分はフローチャートの下半分を占め、密度汎関数部分
(DFT)は上半分を占めている。分子動力学法は零K
でない所謂有限温度でのやや大きな系を扱い、原子(イ
オン)間のポテンシャルには、予め第一原理から求めた
ポテンシャルを用いる。他方、密度汎関数(DFT)
は、小さな系の基底状態を求める時に使うことにした。
そして、当然DFTの部分は温度は零Kとした。ポテン
シャル部分は、電子系の基底状態から求めた。具体的に
は波動関数(KS方程式)を解き、局所密度汎関数を併
用した。
【0101】DFT部分で、縮重のない電子系の基底状
態の全エネルギーを求めて置く。全エネルギーは(波動
関数までさかのぼらずに)電子密度n(r)の汎関数E
[n(r)]である。
【0102】即ち、E[n(r)]は、下記式(120
1)により表される。
【0103】
【数1】 右辺第1項のVz (r)は核のクーロン場などの外場、
第2項は電子間のクーロン相互作用エネルギー、第3項
は、電子の1体近似を想定して、電子間相互作用のない
時の運動エネルギーTs[n(r)]であり、下記式(1
202)により表される。
【0104】
【数2】 第4項は交換相関エネルギーExc[n(r)]で、その
形は局所密度近似による。電子密度n(r)は下記式
(1203)により表される。
【0105】
【数3】 基底状態は、式(1201)のE[n(r)]の変分を
0とおいて得る、下記式(21),(22)に示すKS
方程式(Kohn−sham equation)を解
いて、固有値Eの小さい方からN個採用する。式(12
11)にn(r)が入っているから、n(r)、ψ
(r)、有効ポテンシャルVeff (r)を自己無撞着
(self−consistant)に解かなければな
らない。
【0106】
【数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を掛けて積分することにより、下記式(121
3)を得た。
【0107】
【数5】 Exc[n(r) ]=∫εxc(n(r))n(r)dr …(1213) εxc(n)の形は、Wigner他の式など、いろいろ
提案されている。更に、内殻電子は考えずに価電子だけ
扱う時には、式(1211)のVz (r)として核のク
ーロン場をそのまま使うのでなく、核の位置での特異性
を消したノルム保存擬ポテンシャルでおきかえる。この
ようにならすことによって、フーリエ展開での成分の個
数をおさえることができた。
【0108】全エネルギーとして、本発明者等は厳密を
期するため、原子(イオン)間のクーロン・ポテンシャ
ルも含むものを使った。即ち、下記式(1221)であ
る。
【数6】 ここで{RI}はイオンの核の座標、{αν}は外部の
制約条件(体積Ω、ひずみεLνなど)。式(122
1)の右辺第1項は、DFTの式(1201)から右辺
のTs [n(r)]であり、式(1201)右辺の残り
が式(1221)のUの1部である。その部分は、DF
Tの場合と同様に1体近似の式(1211)からVeff
(r)で表され、やはりVz (r)を擬ポテンシャルで
置き換えて使う。
【0109】ここで新しい考え方を入れた。即ち、Eを
ポテンシャル・エネルギーと見なして、他に、下記式
(1222)に示す仮想的な運動エネルギーを導入す
る。
【0110】
【数7】 従ってラグランジアンは、下記式(1223)により表
される。
【0111】 L=K−E …(1223) 式(1222)中のMIは、本当の原子(イオン)の質
量だがμ、μνは仮想的な質量であり、後述のように目
的に応じて変えた。束縛条件として、下記式(122
4)に示す正規直交条件
【数8】 ∫ψi *(r,t)ψk(r,t)dr=δik …(1224) を課するから、オイラー方程式を作る時に未定乗数Λik
を導入して、Lに下記式(1231)を加えたものの変
分をとる。その結果、下記式(1232),(123
3),(1234)が得られる。
【0112】
【数9】 式(1222)の運動エネルギーの温度Tが対応する。
古典統計力学のエネルギー等分配則によって各自由度が
温度を持っているのだから、Σ1/2 MII 2 との関係
でいえば、仮想的でなく現実の温度である。またμ、μ
νの大小によってψi 、ανを抑制したり増強したりす
ることが出来る。
【0113】例えば、μ<<MIとして高温からT→0と
すれば、RI をとめたままψi が変化して電子系の基底
状態を得る。この時、式(1231)の左辺が0にな
り、右辺は式(1221)とその下の注意により、下記
式(30)となるから、ψi の線形結合を適当に作る
{Λik}が対角線化されてKS方程式(1235)を得
る。
【0114】
【数10】 [∇2−Veff(r)]Ψj(r)+ΣΨk(r) …(1235) つまりDFTでKS方程式を解いたのと同じ結果をsi
mulated annealingで得たことにな
る。ただし、温度Tの状態を作るのに、モンテカルロ法
によらず、仮想的な力学系の運動を追っているから、d
ynamicalsimulated anneali
ng(DSA)と呼ばれる。
【0115】式(1231)の積分をする時、必ずしも
ψi をそのまま変数として扱うのでなく、それを平面波
で展開した展開係数(つまりフーリエ係数)を変数とし
て扱うことができる。この展開係数の個数をむやみにふ
やさないためには、DFTの場合と同様に、擬ポテンシ
ャルで核の位置での特異性を消しておくことが必要であ
る。特に、結晶の周期性を利用できるときには、上記の
{Λik}を対角線化するψi の線形結合をφnkと記す。
kはブリユアン帯域内の波数ベクトル、nは帯域インデ
クスである。
【0116】φnkは、下記式(1241)を満たすはず
である。
【0117】
【数11】 μSnk(r,t)=[∇2−Veff(r)]φnk(r)+λnkφnk(r,t) …(1241) λnkは行列{Λik}の固有値であり、エネルギー準位で
ある。φnk(r,t)を展開して、下記式(1242)
が得られる。
【0118】
【数12】 擬ポテンシャルが局所的であるならば(この条件は一般
に満たされない)、有効ポテンシャルVeffも展開し
て、下記式(1243)が得られる。
【0119】
【数13】 式(1242)と式(1243)を式(1241)に代
入して、下記式(1251)が得られる。
【0120】
【数14】 右辺の最後の頃でK=K′+K″とおくと、Σexp [i
(K+G)r]が全部の項に共通となるから外して、下
記式(1252)を得る。
【0121】
【数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)の左辺を差分になおすと、下記式(12
61)が得られる。
【0122】
【数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を大きくすることに
より、計算量を減少させることが出来る。
【0123】クーロンポテンシャルは、以下のようにな
る。即ち、本発明者は、クーロンポテンシャルを分解し
てゆくと4項に分かれる事を見い出した。特に、従来は
3項しかなかったが、新しく第4項があることを確認し
た。これらの方法を順次説明する。
【0124】本発明者らは、まず厳密な式の出発とし
て、誘電率を加味した基本式から解き始めた。まず基本
式である下記式(1262)から解き始めた。
【0125】
【数18】 導体ε=∞と真空ε=1の場合とでクーロンポテンシャ
ルが異なり、Lは(立方体の)単位結晶の一稜、Σは単
位結晶内でとり、Zi 、ri は、それぞれ第i粒子の荷
電と位置である。これは球内の荷電によって球の内面に
分極が生じることによる。導体でない球の内面に双極子
の層が出来るが、上記式の最終項がちょうどその効果を
打ち消す働きをする。Ewaldの方法は左辺、ε=∞
のものを与えるから、真空内の値を求めるには上記式の
最終項を加えなければいけない。結果だけ掲げる。
【0126】クーロン・ポテンシャルを下記式(127
1)により表わすものとする
【数19】 式(1271)において、Nは単位結晶内の原子数、そ
の中の第1原子の位置と荷電がri ,Zi であり、nは
単位結晶とそれを周期的にずらしたものを指定するベク
トルであって、下記式(1272)に示すように設定し
た。
【0127】 n=nxx+nyy+nzz …(1272) と表される。ここに、Zx ,Zy ,Zz がそれぞれx,
y,z方向の稜のベクトルで、nx ,ny ,nz はそれ
ぞれ(バルク結晶の場合)−∞から+∞にわたる整数で
ある。Σ′の′はn=0の時のj=iを除くことを意味
するものとする。
【0128】ここで本発明者らは下記式(1273)に
示す新しいF関数を導入する。
【0129】
【数20】 上記式(1273)において、G(r,t)はtの周期
関数である。これはフーリエ級数で表現できることを見
い出した。
【0130】G(r,t)は、更に下記(1281)に
示すように変形することが出来る。
【数21】 上記式(1281)において、mは逆格子ベクトルであ
り、下記式(1282)により表わされる。
【0131】
【数22】 m=mx・(1/Lx,0,0)+my・(0,1/Ly,0) +mz・(0,0,1/Lz) …(1282) 指数mx ,my ,mz はそれぞれ−∞から∞にわたる整
数である。m=0の時は、exp […]=1だから、ΣZ
i =0であり、単位結晶内の総荷電が0であれば、m=
0の項は消える。あらかじめm=0の項を除くと、下記
式(1283)となる。
【0132】
【数23】 本発明者等は積分範囲をkで分け、G(r,t)の2つ
の形を使い分けて、下記式(1291)を得る。
【0133】
【数24】 ここで本発明者等は、公式として下記式(1292)を
用いた。
【0134】
【数25】 これにより、始めのクーロンポテンシャルは、下記式
(1293)により表わされる。
【0135】
【数26】 Σ′の右上の′印は、n=0の時、j=iを除くことを
意味する。
【0136】更に、下記式(1301)が成立する。
【0137】
【数27】 この式(1301)はrを含まないから、力には関係な
い。この式(1301)は、更に下記式(1302)に
示すように変形することが出来る。
【0138】
【数28】 以上の結果をまとめると、クーロンポテンシャルは、結
局、下記式(1311)〜(1315)により表わされ
る。
【0139】
【数29】 ここに、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) での減衰
との遅速に逆向きに効くから、両者のバランスがとれる
ような適当なκを選ぶ必要がある。これらは、距離の近
いところから順にクーロン力の寄与を計算した結果であ
り、しかも周囲が導体である場合の結果である。
【0140】周囲が真空である場合によればもう1項が
加わる。これが特に従来省みられなかった部分である。
【0141】一般に分子動力学では、原子間ポテンシャ
ルが最も大切な量である。これには、現状では、単純に
は、2体間でのみ記述するものも結構多い。これは、分
子動力学計算が非常に、時間のかかるもので、その時間
節減のために、ポテンシャルを簡便なものを使っている
わけである。この様に、ポテンシャルを簡便にしてしま
うと、計算自体は速くなるが、実際を反映したものでは
ありえない。
【0142】上に記したように重要な原子間ポテンシャ
ルは、一般には、3体で記述できる。これの意味合いは
以下の通りである。即ち、第一の原子(i)と、第二の
原子(j)との作用に、さらに第三の原子(k)から第
二(j)を通して入ってくる効果がある。これを図示し
たものが図4である。図に有るように、rを粒子間の距
離とし、θを粒子角度とする。これらの諸量の主従の関
係は、図5に示す通りである。特に図5で分かるよう
に、非常に複雑であることが分かる。一般に、ポテンシ
ャルを微分したものが、力になるわけであるが、これで
分かるように計算は、非常に複雑になる。ここでもう一
度、従来の計算手法と本発明の位置つけを図6を用いて
示した。
【0143】本発明者らは、図6にある高階巡回偏微分
法を新しく開発して、これにより、高精度に且つ高速に
演算できる様にした。本発明者等の方法がどれくらい速
いかを以下に実際の比較例を入れて記す。まず、図6に
ある簡便な2体間ポテンシャル計算に関しては、ここで
は、触れない。なぜなら、半導体LSIの物理現象を記
述するには殆ど使用に耐える精度が得られないので、こ
こでは、議論しないことにした。一方、3体問題を入れ
たポテンシャルで、従来の方法として、テーブル参照方
法のプログラムを、本発明者等は作成した。これによる
と、rijと,θi jkと、rjkに関して、表を作成した。実
際には、block dataとして、配列にして格納した。そし
て、この参照項目からずれていた時のことを考え、補完
プログラムをも作成した。補完には、cubic spline法を
用いた。ある値、即ち、rijと,θijkと、rjkとに関し
て、テーブルを引用出来るようなプログラムにした。し
かし、所望の条件(rijと,θijkと、rjk)が、配列内
のどの行列に入るかを探る必要がある。本発明者等は、
出来るだけ高速に成るように気を配りながらプログラム
を作成したが、やはり、数字の大きさによる場合分けを
する部分が必要で、そこにはif then else文を用いた。
この様にして、所望の条件(rijと,θijkと、rjk)が
テーブルのどの枠にはまるかを、見いだすとともに、そ
の中での補完をおも行った。これで、まず、所望の条件
(rijと,θijkと、rjk)でのポテンシャルが求まった
わけであるが、実は欲しいのは、この点における力であ
る。本発明者等は、微小変化させた部分でもう一点につ
いても、ポテンシャルを同様の手続きを用いて求め、両
者の差から、平均変化率の考えを用いて、微分、即ち力
を求めた。
【0144】これで分かるように、非常に多くのif the
n else手続きを用いた。電子計算機では、if then else
手続は比較的苦手で、遅い。事実、本発明者等が作成し
た、参照型プログラムを用いた場合、力を計算する部分
のみについて、詳細について調べた。その結果、160
0個のSi原子からなる単結晶基板で、絶対温度300
Kで、外圧1気圧の場合であるが、参照型プログラムで
は、約350倍の時間を要した。これで分かるように、
本発明者等が提案する、高階巡回偏微分法のほうが、有
利であることがわかる。また参照型プログラムでは、原
子が異なる度に、その都度データ表を作成する必要があ
る。
【0145】次に分子動力学から材料特性値等をどの様
に抽出して、それを汎用シミュレータにどの様に入れて
いるのか」に関しても簡単に説明する。
【0146】通常、LSIプロセスに於いては、種々様
々な評価技術を使っている。赤外吸収による膜の物性値
評価法、さらには、ラマン分光による膜中の応力評価、
X線による膜の構造解析等々がある。またこれら光学的
なものに限らず、膜の誘電率等電気的な評価技術もあ
る。
【0147】本発明者等は、それぞれの評価手法の原理
を、一旦、鋭意、理論的に個々原子の運動にまで分解し
詳細に追跡した。その上で、本発明者等の開発した分子
動力学を用いて、個々原子の動きを追跡し、それぞれの
評価量を初めて理論的に算出し、たとえばスペクトルの
様な形に直接表示することが出来た。これにより、例え
ば、実測値との定量的な比較が出来る等、従来にない大
きな技術の進捗を得たわけである。
【0148】本発明者等が、初めて克服した理論的手段
による原子運動の追跡と、具体的な測定評価量の導出ま
でを以下に記す。幾つかのものがあるが、基本的には、
本発明者が導出した基本概念は以下の通りである。即
ち、着目する原子のモーメントを計算する。そして、こ
れを、着目する時間の窓の中で積算し、しかる後に、フ
ーリエ変換し、固有スペクトルを得るものである。これ
が大きな骨子である。
【0149】次に、Siを例に、経験的なポテンシャル
に関して、本発明者がどの様にして、そのポテンシャル
を力に変換したかを示す。これは、一例であるが、Si
に限らず、イオン結合による効果を考慮しなければ、例
えば酸化膜のポテンシャルも、このSiの問題と同様
に、3体問題として、ここに述べる高階巡回偏微分を用
いることが出来る。本発明者等は、Siに関してどんな
ポテンシャルが良いかも吟味した。
【0150】ここで用いるポテンシャルはTersoffによ
るものである。(Phys.Rev.Lett.56.632 (1986), Phy
s.Rev.B37,6991(1988))。
【0151】今、Tersoffによれば、i番目のSiに関す
る全ポテンシャルは、
【数30】 で記述できる。これは3体で記述しているので、本発明
者等が注意を喚起したいのは、上記式(0001)に於
いて、Vij≠Vjiである。着目するSi原子の位置番号
をiとし、その周辺の他の粒子番号をjとすると、上記
ijは、
【数31】 Vij=fc(rij)[aijR(rij)+bijA(rij)] …(0002) である。ここにrは粒子間の距離である。また、fc
(rij)は、カットオフ関数で、fR(rij)は斥力を
示し、fA(rij)は引力を示す。aijは配位数を考慮
したカットオフ係数、bijも配位数を考慮したカットオ
フ係数である。配位数を表現するパラメータを用いて3
体的表現を行っていることになる。これは、他のもので
も基本形は同じであるので、ここでは、本発明者等が初
めて開発した巡回手法の思想を交えながら以下に示す。
【0152】fR(rij)とfA(rij)は、Morse型の
ポテンシャルを変形したもので、fR(rij)=Aexp(-
λ1r)、fA(rij)=−Bexp(-λ2r)である。
【0153】この内、λ1とλ2は、定数であり、その大
きさは、原子間距離程度の値の逆数である。
【0154】ところで、カットオフ関数fc(rij)と
は、
【数32】 Vij=fc(rij)[aijAexp(−λ1ij)−bijBexp(−λ2ij)] 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Åである。次に、b
ijは実配位数を示すことになるが、ここでも上記カット
オフ関数を使う。その定義は、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である。
【0155】また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) である。
【0156】これらの準備を経ていよいよ本発明者の巡
回部分を詳しく説明する。ここからいよいよ力を計算し
て行くわけである。ポテンシャルの(0002)式を位
置の座標で微分すると力になる。即ち、
【数38】 がそれぞれ、粒子i、jに働く力のベクトルのx成分であ
る。しかし、実際にそれを求めるのは、そう簡単ではな
い。本発明者等の数学理論に立ち戻った工夫と発明が以
下に有るわけである。それが、高階巡回偏微分と言う分
けである。
【0157】(0011)式及び(0012)式で述べ
た式を、具体的に偏微分すれば良い。しかし、そんなに
単純では無いことが直ぐ分かる。
【0158】即ち、
【数39】 また、j及びkに関する偏微分方程式の変形は、同様に
以下の様になる。また、上記は主にxについて記述した
が、yやzについても行う必要がある。ここでは、特
に、上記との対応が分かるように、空白部分は、空白の
ままにして置いた。このような複雑な変形は、本発明者
等が、調べた限りでは見あたらなかった。この様な複雑
な変形を今まで克服しなっかたため、逆に、単純なテー
ブル参照方式に逃げたりしていたわけである。この様な
発明を克服したため、本発明者等は、分子動力学と汎用
流体シミュレータとを合体出来たわけである。
【0159】
【数40】 これらの式の各項は次のように変形できる。
【0160】
【数41】 ∂Vij/∂rij =∂fc(rij)/∂rij・[Aexp(−λ1/rij)−bij Bexp(−λ2/rij)]+ fc(rij)[−λ1 Aexp(−λ1/rij)+λ2ijBexp(-λ2/rij)] =Aexp(−λ1/rij)・[∂fc(rij)/∂rij−λ2c(rij)]− Bexp(−λ2/rij)・[∂fc(rij)/∂rij−λ2c(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 3c(rik)fc(rij)(rij−rik)2 …(0020) ∂ζij/∂cosθijk =fc(rik)exp(λ3 3(rij−rik)3)dg(θijk/dcosθijkc(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/(rijik)・∂Pijk/∂xi +ijk{1/(rik)∂/∂xi[1/(rik)]+ 1/(rij)∂/∂xi[1/(rik)]} =1/(rijik)・(xi−xk+xi−xj)− Pijk{(xij)/(rikij 3)+(xi−xk)/(rikij 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)等から転送されるものが
ある。内容は、境界条件や、領域の大きさ(即ち、x方
向の広がり、y方向の広がり、z方向の広がり)、応力
(即ち、σx、σy、σz)、歪み(即ち、εx、ε
y、εz)、温度などである。また必要に応じて第二の
ものとして補助的な入力bがあり、分子動力学の計算時
間やアンサンブルの指定等がある。
【0161】他方、第一の出力としては、位置ベクトル
r(t)が出てくる。そして、図8に示したファイル
(23)、(33)、(43)、(53)、(63)、
(73)、(83)等へ転送される。叉、膜特性算出プ
ログラムIIIにより、光学的・電気的・結晶学的特性
等の膜特性を算出し、図8のファイル(3)へ転送され
る。
【0162】また、図9に、汎用2・3次元プロセスシ
ミュレータ(図8の(Ia))の入力シーケンスの一例
を示す。また、図10に、同じく汎用2・3次元プロセ
スシミュレータ(図8の(10))のもうひとつの入力
シーケンスを示す。図9に示した内容は、例えば、それ
ぞれ用いるクラスタツールをも指定するとともに、例え
ばその中にある酸化工程の詳細には、ガス量や、昇温シ
ーケンス、さらには、バルブシーケンス等がある。
【0163】他方、図10の入力シーケンスは、図8で
示した未知材料を含む場合である。この図10の場合
は、例えば、強誘電体膜であり、膜をスパッタする工程
を含んでいる。このとき、強誘電体膜の詳細データがな
いので、分子動力学シミュレータIIに、一旦、飛び、
そこで計算を行われる。
【0164】又、本発明者らは、経験的ポテンシャルを
与えた分子動力学シミュレーションを行い、(1)各プロ
セス要素の揺らぎや各要素のずれの重畳現象も診断し、
(2)所望のシーケンス通り実行されているかのチェック
や、(3)変動因子の発見も可能とするシミュレーション
システムを開発した。
【0165】具体的には、単結晶強誘電体膜とSi基板
との位置関係や、金属電極との位置関連を提案するとと
もに、さらに、単結晶強誘電体膜として、実際上どの程
度の品質のものを用意すれば良いかについても、初めて
指針を示すことが出来た。かかる単結晶強誘電体膜構造
を厳密に独自の計算手法により、計算機を用いて再現す
るとともに、特性評価関数を付与し、これに基ずき、印
加条件を考えて基板単結晶との位置関係をどの様にすれ
ば良いかと同時にその欠陥密度の許容範囲の指針を示す
ものである。
【0166】本発明者らは新しい計算化学理論により、
種々元素を含む系における様々な現象を高速に且つ高精
度に計算予測する方法を発明した。これによって、半導
体製造プロセスに使用される様々な材料の特性が高精度
に計算できるようになった。
【0167】従来、種々元素が含まれている結晶や分子
の構造及び特性を予測する手段としては、abinit
io分子動力学法や第一原理分子軌道法が使用されてき
た。abinitio分子動力学法は主に結晶の構造や
電子状態の調査に使用され、分子軌道法は通常の分子
や、結晶の一部を切り出したクラスタを対象として、そ
の構造や電子状態の計算を行うのに使用されてきた。と
ころが、これらの方法では、系の各電子に対してシュレ
テ゛ィンガー方程式を解くため、膨大な計算時間が必要で
あった。そのため、経験的ポテンシャルを与えた分子動
力学法のように、数10万ステップもの計算回数を必要と
する手法の中で、これら第一原理手法を取り込む事は現
実的に不可能であった。一方、分子動力学の分野では、
1991年に、A.K.RappeとW.A.Goddardが"Charge Equilli
bration for Molecular Dynamics Simulations"と
いう論文の中で、多原子分子に対する新しい電荷計算手
法を提案している。本発明者らは、A.R.Rappeらによっ
て提案されている多原子分子における電荷計算手法をさ
らに大きく拡張し、種々の元素や結晶系に適用できるよ
うに大幅に変更し、半導体装置の材料として一般に用い
られるSiやSiO2はもちろん、PbTiO3やPZTといった強誘
電体や、高誘電体にも適用できる材料設計シミュレーシ
ョン手法を構築した。まずはじめに、Rappeらの手法に
ついて以下に述べる。
【0168】Rappeらは、まず、孤立原子のエネルギーE
をイオン化ポテンシャルと電子親和力の関数として考え
た。すなわち、ある原子Aにおいて、電気的に中性な状
態をA0で示す。テイラー展開の2次の項まで考えると、
【数44】 原子が電気的に中性な状態のエネルギーは、(3031)
式のQAに0を代入して得られる。
【0169】 EA(0)=EA0 …(3041) 原子が+1価イオン化した時のエネルギーは、(3031)
のQAに1を代入して得られ、-1価にイオン化した時のエ
ネルギーは、-1を代入して得られる。
【0170】
【数45】 (3042)と(3041)の差はイオン化ポテンシャルIP
になり、(3043)と(3041)の差は電子親和力EAに
なる。
【0171】
【数46】 (3045)式のEAに−が付いているのは、電子を取るの
を正の仕事と考えると、電子を加える事は負の仕事をし
た事になるためである。(3044)と(3045)の差を
とると、電気陰性度χA0となる。
【0172】
【数47】 また、(3044)と(3045)の和をとると、次式が成
り立つ。
【0173】
【数48】 ここで、(3052)式の物理的意味を考える。水素原子
の場合のイオン化時のクーロン相互作用の模式図を図1
2乃至図14に示す。図から分かる様に、IPとEAの差は
電子間クーロン相互作用bになっている。従って、φA軌
道中の2電子間のクーロン斥力をJAA0と表すと、 IP−EA〓JAA 0 …(3053) JAA0をidempotentialと呼ぶことにする。電子を加える
事により、軌道形状が変化するので、Hartree-Fock等の
第一原理計算で計算されたJAAと(3053)式の値は多
少異なってくるが、ここでは(3053)式で近似する。
【0174】Idempotential JAA0と原子直径RA0の関係
を考えてみる。図15に示した簡単な原子構造モデルを
考えると、電子間のクーロン相互作用JAA0は以下の様に
なる。
【0175】
【数49】 ここで、 JAA0はeV単位、RはÅ単位である。(3053)
式の妥当性を確認するため、J値とRとの関係を調べ、表
1に纏めた。表1に示した様に、Jは原子サイズにほぼ逆
比例し、(3054)式から求められた原子の直径は等極
結合距離とおおよそ一致していることが分かった。
【0176】(3051)(3053)式を用いると、孤立
原子のエネルギー(3031)は、次式で書ける。
【0177】
【数50】 多原子分子の場合には、原子間の静電エネルギーも考慮
する必要がある。静電エネルギーは次式で表せる。
【0178】
【数51】 ここで、JABは原子A,Bの中心に1電子電荷がある時のク
ーロン相互作用である。JABはAB間の距離に依存する。
(3055)式を用いると、全静電エネルギーは(305
7)式で表せる。
【0179】
【数52】 原子スケールの化学ポテンシャルχI(Q1,・・・・,QN)は次
式で表せる。
【0180】
【数53】 電荷平衡の条件(3059)から、(N-1)ヶの連立方程式
が得られる。
【0181】 χ1=χ2=………=χN …(3059) すなわち、
【数54】 となる。これに、全電荷の条件(3062)式を含めて、
N個の連立方程式を解く。
【0182】
【数55】 (3061)(3062)をマトリックスで表示すると、以
下の式で表せる。
【0183】
【数56】 (3063)式を解くに当たって、各イオンが取り得る電
荷には、上限下限があるので、電荷がその範囲外なら
ば、その電荷を境界値に固定する。(3061)式を境界
値に固定した電荷と、そうでない電荷とに分離し、(3
063)式を修正してマトリックスを解き、収束させ
る。すると、電荷のばらつきが計算できる。
【0184】この方法の中で問題なのはJABの決定であ
る。原子AとBとの距離RABが大きいときには
【数57】 でよいが、電子雲が重なる距離においては斥力を考える
必要がある。ここでは、1つのスレーター軌道の項で電
子密度を近似してしまう。すなわち、外側の原子価軌道
がns,np,あるいはndである原子において、何れの場合に
も次の型の規格化されたnsスレーター軌道からなるとす
る。
【0185】
【数58】 (3065)式を用いて原子の平均サイズを求めると、
【数59】 となる。下記の関係によって、原子価軌道指数ζAを選
ぶ。
【0186】
【数60】 λは調整用のパラメタで、(3066)式で与えられた平
均的な原子サイズと結晶の共有結合半径RAとの差を計算
するために含まれている。
【0187】スケーリングパラメタλは、図16に示し
たようなアルカリ金属ハライド分子を用いて設定する。
アルカリ金属ハライド分子をMXとし、Mをアルカリ金
属(Na、K、Rb、Cs)、Xをハライド(Cl、B
r、I)とする。それぞれの電荷をQM、QXと表すと、(3
062)式から、 Qtotal=QM+QX=0 ∴QX=−QM …(3068) (3061)式から
【数61】 (3069)式に(3068)式を代入すると、
【数62】 (JMM−JMX−JXM+JXX)QM=χ0 x−χ0 M …(3071) ところで明らかに JMX=JXM …(3072) であるから、
【数63】 (3073)式の変数はλだけである。計算されたQから
(3074)式を用いて双極子モーメントμMXを求め、実
験値と比較してλを決定する。
【0188】 μMX=4.80324QMMX …(3074) RMXは結合距離の実験値である。(3074)式はQがe単
位、RがÅ、μがデバイの単位になっている。双極子モ
ーメントの実験値との比較から、λの最適値として0.49
13が得られている。このλを用いると(3067)式より
原子価軌道指数ζが図17の表の値に求まる。
【0189】水素に関しては、Rappeらは、他の元素と
異なる取り扱いをしている。
【0190】水素の電気陰性度をMullikenとPauling他
の実験上の値と比較すると一致していない。Paulingスケー
ルでは、水素は炭素より電気陽性で、ホ゛ロンよりわずかに
電気陰性であるのに対し、Mullikenスケールでは、水素は炭
素や窒素より電気陰性である。この電気陰性度の問題は
結合している水素の軌道は自由な水素イオンのように広
がる事ができないので、実行的な電子親和力EAが原子の
時の値よりも小さくなるためである。
【0191】従って、EAHを変数にして、水素のχ0HとJ
0HHとを再定義する。
【0192】
【数64】 実験との比較から、最小二乗法によって、 χ0 H=4.528eV J0 HH(0)=13.8904eV …(3076) 実験でなく、HF波動関数の静電ポテンシャルから計算さ
れた電荷との比較から、χ0 H=4.7174eV J0 HH(0)=13.4725eV …(3077) を使用しても良いとしている。
【0193】また、電荷に関する境界値を設定する。
【0194】χA 0,JAA 0は、明らかに、原子価殻が空
っぽ、あるいは満杯に対応する範囲外では妥当でない。
例えば、Li,C,Oの場合には、下図に示したような
範囲に電荷を制限する。この範囲外ではEA(Q)=∞
をとる。
【0195】マトリックス(3063)を解いて求めた電
荷を、図18乃至図20に示した不等式の範囲内にある
か否か判定し、その範囲外である場合には、その電荷を
境界値に固定する。
【0196】以上がRappeらの手法のあらましである。R
appeらは、前述の手法を有機系の分子に適用している。
本発明者らは、Rappeらの提案した原子スケールの化学ポテ
ンシャルという概念だけを採用し、結晶系における個々
原子の電荷計算手法を新たに開発した。以下に本発明者
らによる電荷計算手法を述べる。
【0197】(3057)式右辺第二項は、すべての結合
の静電エネルギーの和を意味する。本発明者らは、分子
から結晶に拡張し、周期境界条件を取り入れようと試み
た。すると、静電エネルギーを計算する際、計算に用い
ている数百〜数千の原子だけでなく、その外側にある無
数の仮想原子からうける静電エネルギーも考慮する必要
があることを発見した。すなわち、静電エネルギーは減
衰が小さく、遠方の原子からの影響も蒸し出来ないこと
に気がついた。そこで、計算に用いているセルの外側に
仮想的に無数のセルが存在すると仮定して、仮想セルか
らの静電エネルギーも考慮するようにした。
【0198】この場合、(3057)式仮想セル中の原子
から受ける静電エネルギーJABは 、計算セルの大きさを考
慮すると、RABが通常10オングストローム以上の大きな
値 となるので、クーロンエネルギーだけと考がえて良
い。従って、JAB=1/RABと一 律に扱う。仮想セルをνで
表すと、(3051)式は次の様に変更される。
【0199】
【数65】 ここで、ν=0は計算セル自体を意味し、ν≠0は仮想セルを
意味する。(3078)式第3項にエワルト゛法を適用する。エワ
ルト゛法を用いると、第3項は次の様に求めることが出来
る。
【0200】
【数66】 Gは逆格子ベクトル、Ωは計算セルの体積、κはエワル
ド計算の収束を調節するパラメタである。更に、(30
82)式右辺第一項はν=0且つA=Bの項は含まず、第二項
はG=0の項は含まない。これらの式から、原子スケールの化
学ポテンシャルχIを求めると、
【数67】 となる。電荷平衡の条件から、N-1個の方程式が得られ
る。
【0201】 χ1=χ2=χ3=……=χN …(3092) (3091)式を代入して、
【数68】 (3093)式で表せるN個の連立方程式を解くことにな
る。
【0202】(3093)式からCQ=-Dのマトリックスを
作ると、次式が得られる。
【0203】
【数69】 本発明者らは、SiO2の様な結晶や非晶質を扱う場合に
は、周期境界条件を考慮しながら(3094)式を解けば
良い事を発見した。(3094)式は電荷Qに対して非線
形となっているため、計算プログラムにはQを収束させ
るためのルーフ゜を持ってくる。
【0204】さらに、本発明者らは、実計算にあたって
のクーロン積分のモデル化にも注意を払った。JABは、
ζAとζBに対するスレーター関数の2原子クーロン積分として
求めた。クーロン積分JAB。はRoothernの方法を用いて
計算した。Roothernらは、クーロン積分JABには、 と
いう標識を使用する。 (3101)式で定義されたハ゜ラメタ
を用いて、クーロン積分は(3102)式で表せる。
【0205】
【数70】 ここで
【数71】 (α,α′,β,β′) =(1+τa)α(1−τa)α′(1+τb)β(1−τb)β′ …(3103) である。さらに、
【数72】 ここで、上式の単位系は原子単位系であり、距離はboh
r、エネルギーはHartreeを使用している。
【0206】Roothernらは1s,2s,2pに対するクーロン積
分の解析式を与えているが、本発明者らは、3s以降の原
子軌道に対するクーロン積分まで計算し、半導体装置で
一般的に用いられるSi及びSiO2が扱えるように、最適な
軌道の選択を行った。
【0207】JABの値は電荷計算に大きく影響する。そ
のためJを算出する際に用いるスレーター軌道や原子価
軌道指数ζの設定は重要である。そこで本発明者らは、
Si-OやSi-Siを扱う場合のJの値について調べた。
【0208】Rappeらは「分極電荷がns軌道にあると仮
定」しているが、種々元素に対してどの軌道を考慮する
のが最適であるかという点は不明である。 Siの主量子
数nが3であることから、Siの電荷が3sにあるとモデル化
した場合と、2sにあると仮定した場合の計算を行い、比
較した。このとき、n=3におけるJの解析式はRoothernの
方法に従い、本発明者らが、別途用意した。例えば、Si
-Siに対しては、次式を用いた。
【0209】
【数73】 図21はSi-Si間のクーロン相互作用Jの計算結果であ
る。この場合には、2sと3sの何れを使用してもJの値に
殆ど変化が無いことが分かる。ところがSi-O間のクーロ
ン相互作用Jは2sと3sを使用した場合で大きく異なって
いる。図22にOの電荷が2sにあり、Siが2sあるいは3s
にあるとモデル化した場合のJABの計算結果を示した。S
iの電荷を3sだけにおくと、3s軌道が非常に広がってい
る影響を受けて、Jの値がr=1.6Å近辺で極小値を持って
しまう。Si-Oの結合距離は丁度この当たりにあるので、
Jが非常に小さく見積もられ、結果、SiとOの何れの電荷
も極めて過小に評価されてしまう事が分かった。Si-Hの
場合も同様であった。従って、Siの電荷のモデルとして
は、敢えて3s軌道で仮定すると、Jの見積もりが全く精
度を欠いてしまう。このことから、スレーター軌道を用
いるのは、電荷の広がりを考慮するためだけとし、モデ
ルとしては、1sあるいは2sのスレーター軌道を用いるの
が最適である事が分かった。これは、Si以降の原子番号
の元素に対しても同様であり、3s以降の広がった軌道を
モデルとすると、原子間距離が近い場合のクーロン力の
計算精度が非常に悪くなることが分かった。
【0210】水素に対する補正にも、本発明者らは、プ
ログラム上の注意を払った。 (3031)式から、
【数74】 水素補正した時の原子スケールの化学ポテンシャルχHl
を求める。まず、(3057)式に対応する多原子分子の
エネルギーの式は、全原子数N個の内、水素M個とする
と、
【数75】 ここで、(3114)式の右辺第一項、第二項、第五項の
Σ′の記号は、それぞれの和を求めるに当たってHは除
いている事を意味している。原子スケールの化学ポテン
シャルχHlは、
【数76】 (3121)式右辺第四項は次の様に解きほぐす事が出来
る。
【0211】
【数77】 (3122)式より、(27)式の第四項と第五項との和は、
Hl以外の全原子との間の原子間クーロン相互作用を意味
する事が分かる。従って、化学ポテンシャルは、
【数78】 (3123)式を用い、水素補正時のマトリックス(30
63)に変更を加える事になる。今、χ1が水素以外の原
子であるとすると、(3059)式に(3123)式を代入
し、
【数79】 マトリックス(3063)と比較すると、水素のJii0を
【数80】 に変更する事になる。
【0212】次に原子間クーロンポテンシャルJAHの水
素補正を考える。水素の場合のスレーター軌道関数は、(30
75)式より、
【数81】 水素が含まれる場合には、電荷を求めた後に(3126)
式を用いてJAHを修正し、再度電荷を求める。そして電
荷の値が収束するまで、これを反復する。この際、I=1
の行が水素の場合、水素でない行と入れ替えてマトリッ
クスを作成するように、本発明者らはアルゴリズムを作
成した。
【0213】電荷境界条件の扱いにも、本発明者らはマ
トリックスの具体的な変更方法について注意を払った。
【0214】(3061)式を、境界値に固定した電荷QB
と、そうでない電荷QCに分けて考える。
【0215】
【数82】 (3062)式を(3131)式に変更すると、マトリック
スの要素は、のとき、
【数83】 上式の様にマトリックスを変更して、電荷を境界条件の
範囲内に収める。この時、I=1の元素が境界値に設定さ
れる場合、マトリックスの行を入れ替え、境界値設定さ
れていない元素に対する行がI=1になるように行の選択
(ピボット選択)をするようにした。
【0216】本発明者らは、結晶系及び非晶質系に新た
に開発した上記電荷計算法を用い、これを分子動力学法
に追加して、個々原子の動きを逐次計算し、双極子モー
メントを時々刻々計算して記録し、これをフーリエ変換
してIRスペクトルの算出を行った。図23は、以上の手
続きをまとめたものである。
【0217】一方、16G以降のDRAMに使用される
ゲート酸化膜やNAND EEPROMのトンネル酸化
膜においては、高信頼化が従来にも増して重要な課題と
なっている。特に、EEPROMにおいては、高電界下
でトンネル電流を流しながら絶縁性を維持させるという
過酷な条件下で酸化膜を使用するため、絶縁材料として
の限界を原子・電子レベルで探索することが重要となっ
ている。
【0218】このような状況から、Si/SiO2界面
構造やSiO2ネットワーク構造の原子レベルでの探索
実験や、SiO2のトラップ機構の実験、及び理論計算
による予測等が種々の研究機関で進められている。しか
し、本発明のように、ネットワーク構造とIRフルスペ
クトルとの関連を原子レベルで関連づけ、利用する手法
は、従来にはなかった。僅かに、IRスペクトルのSi
−Oの非対称伸縮ピークと酸化界面の応力とを関連づけ
る実験結果が報告さているだけであった。
【0219】また、実験結果の解釈としても、「中心力
モデル」を用い、現実のネットワークを非常に簡略化し
たSi−O−Si分子構造の振動解析からの議論しかな
されておらず、非晶質特有のSi−O結合長の分布、S
i−O−Si、O−Si−O結合角度の分布等を全て考
慮した議論は、従来全く行われていなかった。本発明に
おいては、非晶質特有のこれら様々な分布を考慮するこ
とにより、中心力モデルとは異なった視点で、電気的特
性と構造との相関が得られたのである。
【0220】ここで、本発明者らが克服した点について
付言する。本発明者らは、古典的分子動力学法を使用す
るにあたり、そのポテンシャルを吟味した上、さらに、
改良を加えた。本発明者らは、最初、常行氏らによって
開発されたTTAMポテンシャルの使用を試みたが、S
i−O非対称伸縮振動が非常に低波数側に表されるとい
う欠点を見い出した。そこで、van Beest等に
よって開発されたBSKポテンシャルを使用した。彼ら
のポテンシャルは、TTAMポテンシャルよりも高波数
側にピークが表れたが、尚、実測値よりも低い波数に伸
縮ピークが表された。
【0221】本発明者らは、このピーク値を実測に近づ
けるため、これら既存ポテンシャルで使用されている、
電荷一定という条件を取り除き、Si−O−Si角等が
大変位した場合の電荷の分布を、古典的分子動力学に取
り込んだ。この電荷分布の計算は、理想的には分子軌道
法を用いて算出できるが、現実問題としては、数百原子
もの構造を一度に取り扱って電荷分布計算を行うことは
現状では不可能である。スーパーコンピュータを用いて
も現実的には不可能である。
【0222】そこで、本発明者らは、A.R.Rapp
eとW.A.Goddard IIIによって提案されて
いる手法を用いて電荷分布計算を行った。引用論文は
“Charge Equillibration fo
r Molecular Dynamics Simu
lation”,J.Phys.Chem.,95,1
991,p.3358−3363である。本発明者ら
は、彼らの論文を参考に数式を展開し、論文中の誤りを
訂正の上、必要なパラメータ算出を独自に行って、彼ら
の手法をSiO2系に適用した。Rappeらの手法の
概略について以下に述べる。
【0223】原子iの個別の電荷によるエネルギーEi
(Q)を、テーラー展開の2次まで考えて、下記式(1
651)で近似する。
【0224】
【数84】 ここで、Xi 0 は電気陰性度であり、下記式(165
2)で表される。
【0225】 Xi 0=1/2(IP+EA) …(1652) Jiiは、原子Iの軌道上の電子間のクーロン斥力であ
り、下記式(1661)で定義される。
【0226】 Jii=IP−EA …(1661) ここでIPはイオン化ポテンシャル、EAは電気親和力
を示す。N個の原子で構成される分子の全エネルギーE
(Q1 ,…,QN )は、個別の原子のエネルギーと原子
間の相互作用静電エネルギーJijQi Qj の和である。
ただし、Jijは原子ijの中心にある単位電荷のクーロン
相互作用である。従って、下記式(1662)が成立す
る。
【0227】
【数85】 原子スケールの化学ポテンシャルは、下記式(166
3)で与えられる。
【0228】
【数86】 この式(1663)を、電荷平衡の条件、X1 =…=X
N に代入すると、下記式(1671)が得られる。
【0229】
【数87】 すなわち、下記式(1672)が得られる。
【0230】
【数88】 さらに、下記式(1673)に示す全電荷の条件を考慮
する。
【0231】
【数89】 ijは既知であって、およそ下記式(1674)であ
る。
【0232】
【数90】 ここでRはI=jのときは原子の大きさ、I≠jのとき
は原子間距離である。従って、行列Cとベクトルdを導
入すると、Qに関する下記式(1681)に示す連立方
程式が成立する。すなわち、
【数91】 と定義すると、 CQ=−d …(1682) となる。論文ではC1i=Qi 、Cd =dとあったが、本
発明者らは、これらを上述の様に修正して使用した。さ
らに、Si−O−H等の原子の大きさ等、必要なパラメ
ータを独自に算出し、連立方程式を解いて電荷分布を求
めた。これによって、ピーク波数が格段に実測と一致す
るようになった。
【0233】分子動力学シミュレータIIで算出された
「戸籍簿」、つまり単位時間毎の結合長(位置ベクトル
r(t))と結合角を示す膨大なデータの一例を、図2
4に示す。このデータは、膜特性算出ユニットIIIに
入力され、ここで基本膜特性(光学的、電気的、結晶学
的特性)が算出される。
【0234】又、プロセス変動算出ユニットIについて
は、その2次元汎用システムの一例が、特開平3−16
4904号に記載されている。
【0235】3次元のプロセスシミュレータは、今回本
件発明者が構築に成功した。以下これを、図8をもちい
ながら説明する。まず技術とシミュレータの構造の概観
をする。本願は、原子レベルに主体を置いているが、他
のプロセス変動算出ユニットIとの新しい合体方法をも
提案しており、より効果的かつ正確で、従来にない新し
い機能をも提供できることをも示している。合体する相
手のプロセス変動算出ユニット(10)は、例えば3次
元プロセスシミュレータIcであり、いわゆる流体モデ
ルに属するものである。とくに、本発明者らが新しく構
築した3次元プロセスシミュレータIcの内容を以下に
記す。要点は2つあり,まず、本発明者らが構築した
「3次元プロセスシミュレータ」は、原子レベルシミュ
レータとの接続が可能であること。次に、この「3次元
プロセスシミュレータ」の部分には、本発明者等が新し
く克服した化学モデルや、数学式を駆使している。これ
らにより、初めて、原子レベルとの合体が可能になった
のである。原子レベルシミュレータから出力されるもの
として例えば、粘弾性係数、応力定数、拡散定数、Ci
j、等があり、これらが、流体的汎用3次元プロセスシ
ミュレータに転送されるものである。図8では、(2
3)、(33)、(43)、(53)、(63)、(7
3)、(83)である。例えば、酸化剤等の拡散定数
は、(23)に転送されるわけである。転送手続きの方
法等は後述する。又、これら粘弾性係数、応力定数、拡
散定数、Cij、等は、逆に流体的汎用3次元プロセスシ
ミュレータから原子レベルシミュレータに転送され、互
いに補完することが可能となっている。つまり、共に3
次元で計算を行ない3次元の結果を互いに交換すること
によって、より精密な予測が実現する。
【0236】汎用3次元プロセスシミュレータ(10)
側で、入力データに従って、数ミクロンレベルの広領域
の工程誘起応力の計算を、(2)や(3)で、行い、こ
の結果を、ファイル(22)やファイル(32)を介し
て、一旦分子動力学シミュレータIIに転送し、あるい
は、再び汎用シミュレータ(10)に戻し、これを繰り
返すことにより、より正確でより多くの機能を利用する
こともできる。計算手続きとしては、例えばある時刻
で、プロセス変動算出ユニットで求めたある程度広領域
の工程誘起応力計算の結果を、一旦分子動力学シミュレ
ータIIに転送し、ここで、物性定数の修正量を計算し
ながら、改めて、原子の動きをみ、再び、修正された物
性定数を汎用シミュレータに送り込み、汎用シミュレー
タで広領域の工程誘起応力等の結果を、進める。これを
繰り返すことにより、もし、十分大きな応力が局所的に
生成されれば、分子動力学シミュレータII側で原子の
動きをみ、塑性変形の可否を見ることもできる。
【0237】勿論上記の汎用プロセスシミュレータ部分
での不純物再分布計算や、応力計算の時には、図8に記
した、点欠陥挙動計算部分(9)が呼ばれている。点欠
陥としてはSi基板中の空格子と格子間Siを扱ってい
る。不純物の拡散にはこれら点欠陥の挙動を正しく取り
入れないと、不純物の再分布は正確には計算されない。
イオン注入計算部分(4)では、もちろん、注入損傷に
よる多量の空格子と格子間Siが生成される。これら
は、熱工程中に再分布するが、その時、不純物とも相互
作用を持ち、不純物分布の計算に反映させなければなら
ない。またシリサイド・形状・応力計算部分(7)で
も、勿論、Si基板中の空格子と格子間Siの挙動に十
分配慮しなければならない。
【0238】ところで、図8のその他(8)について
も、少し触れる。LSI素子試作工程で、酸化・拡散・
イオン注入等は、既に多くのデータを基に、このプロセ
ス変動算出ユニット(10)には、各必要物性定数が格
納されている。例えば、Siの酸化に関しては、後に掲
げる式(1130)や(1131)のPB1やPB2等
がそれである。また、SiO2膜中の酸化剤の拡散定数
などは、後に掲げる式(1531)のD’などがそれで
ある。一般には、この様に、概ねデータがあるものであ
る。
【0239】しかし、LSI開発の中では、必ずしもデ
ータ等が揃っているとは限らない。たとえば、今、強誘
電体膜を用いたメモリーセルの試作開発を考える。CD
V或いはスパッタ工程を用いて、PbTiO3を形成す
る場合を考える。現状の汎用プロセスシミュレータで
は、多くはまだPbTiO3膜中での種々の不純物の拡
散係数も、さほど分かっていないし、また、その応力定
数もさほど分かっていない。また、界面での不純物の偏
析定数もさほど分かっていない。この様な場合、まず
(8)の工程に入り、すぐさま、(84)に示す様に、
分子動力学シミュレータIIに一旦行き、そこで、必要
な元素のポテンシャル等を算出し、それを基に、分子動
力学部分が動き、そして、物性定数を算出する。このあ
とは、既に述べてきた部分と非常に良く似ているが、
(83)からこれらの物性定数がファイルに格納され
る。そして、プロセス変動算出ユニットと同じようにし
て、計算部分が進められる。
【0240】ここで、図9を参照して、具体的な処理に
即した説明を行う。ここでは、酸化工程を取り上げる。
即ち、ウエハーの投入を行い、これを基板処理工程に移
し、しかる後に、酸化工程に入る。
【0241】この酸化工程では、例えば、高純度酸素ガ
スの量を例えば、Poxとし、また、これの分散をσtoxと
して構成する。叉、昇温シーケンスはTempとし、また
これの分散をσtemp、バルプシーケンスをtox、σtox等
の様に構成する。特に、バルプシーケンスは、ここで
は、例えばパルプを開けている時問をさしているわけで
あるが、詳細は実施例で述べたようなものである。
【0242】たとえば、図8では、(2)の酸化・形状
・応力計算部分をアクセスし、ここに分圧Poxや、温度
Temp等が入力される。そして、酸化成長計算や、同時
に応力計算も進められる。この時、(21)を用いて、
分散・変動計算も行う。ここではσpoxやσte即が入力
して、中心値からの変分を効率よく計算するものであ
る。
【0243】これらの計算結果は、例えば、酸化膜の成
長量や、酸化膜とSiの界面の形状や、さらには、不純
物量も計算される。さらに、本発明者によれば、応力も
求められるわけである。また分散・変動からも同じよう
に Pox=Pox+σpoxや、Temp=Temp十σtemp等の部
分についても、酸化膜の成長量や、酸化膜とSiの界面
の形状や、さらには、不純物量が計算される。そして、
これらが、一旦メモリまたはファイル(22)にたくわ
えられ、そして、分子動力学部分に送られる。分子動力
学部分では信号aとしてこれらが入力され、位置ベクト
ルr(t)が算出される。
【0244】図10は、図9と同様の行程図を示してい
るが、ここでは強誘電体膜の形成に関するものである。
【0245】ところで、ここで本発明での診断部分に関
する説明も行っておく。図11に、本実施例による2進
木の一部を示す。ここでは酸化工程に着目した2進木を
示している。図中にも記したように、ここではまず、in
-situモ二夕から送付されるSiO2膜厚と、計算から予
測されるSiO2膜厚を比較する。そして、膜厚が許容
範囲に収まっている時を取り上げ議論している。そうす
ると、次の手続きとして、inlituモニタから送付される
1Rスベクトルカープと、計算から予測される1Rカー
ブとを比較する。そして、図に示した様に、まず、スペ
クトルを比較して、少しずれていた場合で、しかも、ス
ペクトルが高波数側にずれた場合を取り上げる。そうす
ると、次の手続きとして、SiO2中に残留応力が大き
く誘起していることが考えられる。それで、この場合、
ゆっくり降温して、残留応力を軽減させることが出来
る。それで、図9に記載したアニール部分を修正して、
図8の(2)に入力される。他方、スペクトルを比較
し、許容範囲に入っていた場合は、当初のプロセスシー
ケンスのまま続行される。又、膜厚が許容範囲に収まっ
ていない場合には、修正不能として、不良処理を行な
う。
【0246】図8に関して、概ねデータの流れがこれで
分かったと思う。では、もう少し、図8の出力に関して
言及したい。即ち、実行結果は、一旦、ファイル(3)
に蓄えられる。ここでは、各工程で求めた、工程に対応
するスペクトルや、その他の光学的・電気的特性結果が
格納される。これらは、あとで、詳細を述べるが、実測
データとしてモニタをしている量と対比できる様になっ
ている。
【0247】一方、本発明者らの1人は、以前、プロセ
ス変動の変分計算を用いたプロセスシミュレータシステ
ムを提案している。即ち、特開昭63−174331号
および特開平3−164904号である。これらの提案
では、流体的な取扱いにより、酸化/拡散、イオン注
入、化学的堆積(CVD)、エッチング、リソグラフィ
等の各工程を取り扱うことができる。例えば、流体モデ
ルによる拡散工程では、不純物総量と電気的活性化不純
物量とを区別しながら、電気効果を取り入れつつ、不純
物の再分布を求めている。
【0248】特開昭63−174331号および特開平
3−164904号に示すシミュレーションシステムの
概略を、図25〜図49を参照して説明する。これらの
提案では、プロセス変動については次のように扱ってい
る。例えば、酸化、イオン注入、リソグラフィの工程を
例にとると、平均の温度(T)とその許容バラつき(δ
T)、平均時間(t)とその許容バラつき(δt)、平
均圧力(P)とその許容バラつき(δP)、更に平均マ
スク寸法(L)とその許容バラつき(δL)、平均ドー
ズ量(D)とその許容バラつき(δD)などが入力値で
ある。
【0249】プロセスシミュレータには、図25に示す
ように、イオン注入、酸化/拡散、リソグラフィ等の各
サブルーチンに分散変動(deviation)を計算
する部分が付設されている。計算手順としては、これら
の中心値について計算を進めながらその変分を計算して
行く方法を採用している。
【0250】この変分計算について若干説明する。まず
拡散工程における温度変分について説明すると、拡散工
程では、拡散係数をD、不純物分布による内部電界を
ε、電気的活性化不純物数をNとして、Fickの第2
法則に基づく方程式を解いて、不純物濃度分布Cを求め
ている。図26を用いて変分計算の手順を示す。ここ
で、図中の左側を正常プロセス(Normal pro
cess)とする。
【0251】非線形拡散方程式は、下記式(0071)
により表わすことが出来る。
【0252】
【数92】 ここでk、k+1は時刻であり、時刻k+1のときの不
純物濃度Ck+1 が未知数である。右辺のBは、下記式
(1072)の第2項に由来するものである。
【0253】 ∂C/∂t→(Ck+1−Ck)/t …(1072) ここでは陰解法を用いているので、上式に現れる係数行
列はすべてDや平均温度T等を用いて表現されている。
正常プロセスの場合は、所定の拡散時間tn+1 番目まで
を適当に分割し、t1 ,t2 ,…,tn+1 まで順次解を
求めて行く。ここではそれぞれの時刻、t1 ,t2 ,
…,tn+1 においてNewton法を用いてこれらを線
形化して計算をすすめている。
【0254】次に、変動プロセス(Process d
eviated)について説明する。ここでは拡散温度
に変動δTが発生した場合を示している。上式の変分を
とり、時刻kにおける濃度の変化量δCK がわかってい
ると、時刻k+1の時の変分δCk+1 は、下記式(10
81)により表される。
【0255】
【数93】 すなわち時刻k+1において正常プロセスの収束時の値
D,ψ,N,C等を保存しておき、これらを用いてF′
を作成し、またその温度Tに対する変化量を用意してお
き、δCk+1 を求める。ここでF′を作成するための行
列要素の一部を図26に示す。拡散係数の導出には、基
本的には空孔モデルを用いている。この場合、温度変化
によってフェルミレベルが変化し、これによって荷電空
孔の濃度が変化し、拡散係数が変化する。従って、まず
その基になる真性キャリア濃度ni に対してはその変分
は、下記式(1082)に示すようになる。
【0256】
【数94】 δni=ni(1.5/T+0.605/kT2)δT …(1082) また、不純物の内部電界の温度による変化も、これらを
用いて表現できる。更に、酸化性雰囲気における拡散係
数の増大現象(Oxidation Enhanced
Diffusion Effect)の変分に対して
も考慮した。酸化工程では線形酸化速度定数(Line
ar rate constant)及び放物線酸化速
度定数(paraboric rate consta
nt)も変化するのでこれも考慮した。これらを用いて
先のF′を作成する。
【0257】ここで、ある典型的な拡散工程を取り上
げ、その1つの拡散時間における収束の様子を図27に
示す。図27の横軸は反復回数を示し、縦軸は残差を示
している。図27では、正常プロセスにおける収束状況
が示されている。正常プロセスでは、3回のNewto
n loopで収束していることがわかる。この収束時
のD,ψ,N等を用いて先に示した変分行列F′を作成
し、温度変分を計算する。図27に示すように、約50
回で収束し、全体の1/10の計算量で済んでいること
がわかる。
【0258】また、図28、図29、図30は、変分の
精度を示す。図28は、正常プロセスとして1000℃
の時の不純物分布を2次元格子にわたって表現したもの
である。図29は、変動プロセスであって、δT=1℃
の場合で図28をもとに計算したものである。他方、図
30は、1001℃を正常プロセスとして打ち出したも
のである。
【0259】これらの図から、図30の値=図28の値
+図29の値なる関係があることがわかる。また、改め
て図30の値をプロットしてみると、図31に示すよう
になり、高温にズレると高濃度側では不純物濃度が下が
る傾向にあり、低濃度側では増大する傾向にあることが
わかる。
【0260】半導体素子を作成する場合、多くのプロセ
スがあり、各プロセスにおいては複数個のパラメータが
ある。例えば、CMOS工程では主なものだけでも、
(i)N−well領域へのイオン注入、(ii)P−w
ell領域へのイオン注入、(iii)LOCOS工程のた
めのSi34のパターニング、(iv)LOCOS工程、
(v)n+ イオン注入、(iv)p+ イオン注入、(vii)
アニール工程などからなる。プロセスの揺らぎを考える
場合、これらの数多くのパラメータの変動とその度数
(或は頻度)を取り扱う必要がある。
【0261】特開昭63−174331号および特開平
3−164904号では、変分の度数と実行の組み合わ
せについて考えている。今、簡単のためにα,β,γの
3つの工程を考えてみる。それぞれにおいてプロセスパ
ラメータの中心値をC(Center)、正に変位した
値をU(Upper)、そして負に変位した値をL(L
ower)と名づけ、それらの度数を、特開昭63−1
74331号および特開平3−164904号から、図
32に示す様に定義してみた。
【0262】工程αとβとを考えた場合、組み合わせ
は、UUやLLなどの計算実行を無視すれば、CU,U
C,CC,LC,CLの5個の組み合わせが考えられ
る。さらにγのプロセスを引き続き経るときは、図中に
示した7つの組み合わせが考えられる。しかし、特開昭
63−174331号および特開平3−164904号
では、これらを全て取り扱うのではなく、UとLの対称
性を考慮し、Lの場合はUから計算するようにしてい
る。
【0263】このシステムの中心部である推論エンジン
の探索の木の一部を、図33に示す。ここでは、酸化・
拡散工程について検討している。図中の“予”は、統計
解析システムによる予測値を示し“実”は、実測値を示
す。
【0264】このとき用いる許容値は、酸化・拡散・エ
ッチング等全工程を経た値である。このときε値は、2
次元プロセス統計解析システムを実行して得られるもの
である。“実”が“予”±εの外に出た場合は、図33
に示すように、酸化工程に異常があった旨記述する。こ
の場合、温度又は時間がどれだけズレていたかを必要に
応じて、統計解析システムのロードモジュールを実行さ
せる。上記検定は、酸化・拡散工程を中心に説明してい
るが、リソグラフィ工程についても行う。リソグラフィ
工程でのゆらぎはSiO2膜の形状の横ズレ、接合形状
の横ズレなどに相当する。
【0265】推論エンジンには、Common Lis
pを用いているが、これは数値処理のみならず文章(L
ist)の処理、形状の処理に適していること、さらに
システムの保守が極めて簡単であるためである。本シス
テムは、推論エンジンから実測値のデータを読むだけで
はなく、統計解析システムのロードモジュールをも実行
している。
【0266】特開昭63−174331号および特開平
3−164904号に示すシミュレーションシステムに
よると、各工程における正常反復計算の収束時の行列の
係数を保存しておき、且つこれを用いて一次変分行列を
作成し、前記個々のプロセスパラメータの変動、及びこ
れらの重畳効果をも考慮することにより、従来の約1/
100の計算量で変動部分が処理できるようになった。
また、たとえばCommon Lispによる推論エン
ジンと上記システムを連動させることにより、プロセス
マージンさらにはプロセス異常の診断さらに最適化まで
出来ることがわかった。
【0267】特開昭63−174331号および特開平
3−164904号では、基本拡散方程式は、各不純物
に対して下記式(1111)に示すように立てている。
【0268】
【数95】 ここで、x,zは2次元空間の水平方向と鉛直(下向
き)方向の座標であり、tは時間を示し、Cは不純物濃
度、Nは電気的活性化不純物量である。またDは濃度依
存性と増速拡散効果を含む拡散係数であり、下記式(1
112)により表わすことが出来る。
【0269】
【数96】 式中、Di は濃度に依存しない基本拡散係数であり、そ
の温度依存性は、下記式(1113)で表現される。
【0270】 Di=BB・exp(−BA/kT) …(1113) BBとBAは不純物の種類による定数、kはボルツマン
定数である。FEはniを真性キャリア濃度として、下
記式(1121)で与えられる。
【0271】
【数97】 またnは1/2[(CD −CA )+(CD −CA2 +4
i 2 ]である。さらにni は、
【数98】 ni=3.87×1016×T1.5×exp(−0.605/kT) …(1122) とした。
【0272】酸化膜成長速度は、放物性酸化速度定数R
1と線形酸化速度定数R2とを用いてそれぞれ下記式
(1123)、(1124)に示すように表現した。
【0273】 R1=PB3・exp(PB1/kT)…(1123) R2=PB4・exp(PB2/kT)…(1124) 上記放物性酸化速度定数をあらためてB、B/Aと記す
と、水平方向の座標xの点で、酸化膜厚をZo 、酸化時
間をt、τを定数として、下記式(1125)に示す関
係がある。
【0274】 これらを用意しておき上記各係数の温度δTに対する変
分を求めてみる。
【0275】δni =ni ・(1.5/T+0.605
/kT2 )δTとしてnの変分は、下記式(1131)
に示すようになる。
【0276】
【数99】 また、下記式(1132)が成立する。
【0277】
【数100】 成長係数の変分は、下記式(1133),(1134)
により表わされる。
【0278】
【数101】 B=R1,A=R1/R2に書き改めると、下記式(1
141),(1142)に示すようになる。
【0279】
【数102】 そして、t+τの変分は、下記式(1143)のように
なる。
【0280】
【数103】 以上説明したシミュレーションシステムは、当時として
は画期的なものであった。しかし、最近のクラスタ化ツ
ールにおける工程制御やモニタ技術には、もはや使えな
くなっている。
【0281】そこで本発明者等は、更に検討を重ね、本
件発明を完成するに至った。ここでは、ペロブスカイト
単結晶を用いた電界効果型MIS素子の強誘電体膜の最
適設計を行い、これに従い、実際に素子を試作し評価し
発明完成を確認した。
【0282】即ち、請求項11によって提供されるの
は、半導体単結晶基板上に電界効果型MIS素子を形成
するにあたり、該MIS素子の強誘電体膜を単結晶と
し、しかも、該半導体基板と該単結晶強誘電体膜を含む
系に於いて、対象特性を表現するabinitio分子
動力学理論に基ずいた評価関数を導入し、これにより所
望の主動作条件下での系の最適解に従った単結晶強誘電
体膜および基板の設計配置とすることを特徴とする電界
効果型MIS半導体装置の製造方法である。
【0283】叉、請求項12によって提供されるのは、
半導体単結晶基板上に、電界効果型MIS素子を形成す
るにあたり、該MIS素子の強誘電体膜を単結晶とし、
しかも、該半導体基板と該単結晶強誘電体膜を含む系に
於いて、その系の対象特性を表現するabinitio
分子動力学理論に基ずいた評価関数として系全体の自由
エネルギを取り上げ、主動作条件下で自由エネルギが最
小になるように前記単結晶強誘電体膜を配置することを
特徴とする電界効果型MIS半導体装置の製造方法であ
る。
【0284】更に、請求項13によって提供されるの
は、Si半導体単結晶基板上に、電界効果型MIS素子
を形成するにあたり、該MIS素子のゲート絶縁酸化膜
を単結晶とし、しかも、該半導体基板と該単結晶強誘電
体膜を含む系に於いて、その系の対象特性を表現するa
binitio分子動力学理論に基ずいた評価関数とし
て系全体の自由エネルギを取り上げ、この自由エネルギ
が、主動作条件下で、最小になるように前記単結晶強誘
電体膜を配置することを特徴とする電界効果型SiMI
S半導体装置の製造方法である。
【0285】更に、請求項14によって提供されるの
は、Si半導体単結晶基板上に、電界効果型MIS素子
を形成するにあたり、該MIS素子のゲート絶縁酸化膜
を単結晶とし、しかも、該半導体基板と該単結晶強誘電
体膜を含む系に於いて、その系の対象特性を表現するa
binitio分子動力学理論に基ずいた評価関数とし
て系全体の自由エネルギを取り上げ、この自由エネルギ
が、主動作条件下で最小になるように前記単結晶強誘電
体膜を配置するとともに、その単結晶強誘電体膜の酸素
欠陥濃度を0.01%以下に抑えたことを特徴とする電
界効果型SiMIS半導体装置の製造方法である。
【0286】更に、請求項15によって提供されるの
は、Si半導体単結晶基板上に、電界効果型MIS素子
を形成するにあたり、該MIS素子のゲート絶縁酸化膜
を単結晶とし、しかも、該半導体基板と該単結晶強誘電
体膜を含む系に於いて、その系の対象特性を表現するa
binitio分子動力学理論に基ずいた評価関数とし
て系全体の自由エネルギを取り上げ、この系の自由エネ
ルギが、主動作条件下で最小になるように前記単結晶強
誘電体膜を配置するとともに、特に局所的なエネルギの
凸凹が偏差値3σ以内に納め、その単結晶強誘電体膜の
酸素欠陥濃度を0.01%以下に抑えたことを特徴とす
る電界効果型SiMIS半導体装置の製造方法である。
【0287】以下にこれらの最適化手順ならびに、試作
結果を記す。
【0288】ここでは、単結晶強誘電体膜として、例え
ばPZT系ペロブスカイトを例に取り、またSi表面を
例えば(100)面とした。この基板面は例えば金属電
極としてFCC構造のCuでも良い。一番重要な問題
は、単結晶強誘電体膜PZT系ペロブスカイトとSi
(100)面とをどのような結晶位置関係に配置すれば
良いかという問題に置き換えられる。ここでは、本発明
の骨子である評価関数として系の全エネルギを取り上げ
た。しかも、この系としては、Si/PZT界面を取り
上げた。本発明者等は厳密な評価関数としてのエネルギ
表現式を予め作成しておき、これを用いることにより、
印加条件下でのPZTとしてのペロブスカイトと、Si
(100)の位置関係は以下の様にするのが最適である
ことを見い出した。即ち、ペロブスカイトのC4v構造
表現のa1軸とa2軸のなす角の中線を、Si(10
0)面上の[110]軸方向に合わせ、且つc軸を傾
け、第1Siと第3Si位置をSi(100)面上のS
iに一致させることが最適値であることを見い出した。
またこの時、PZT中には0.01%以内の酸素欠損が
許容範囲であることも同時に示した。この様に作成した
本発明の提案する概念図を図34に示した。
【0289】図34に付いて詳細に説明する。図34a
の上半分は、Si側の単一格子の平面図を示しており、
下半分はPZT面のスケッチである。また、図34c
は、[110]Si方向に対して、ペロブスカイトを結
ぶ線を並行に配置させるのが最も良いと指摘している。
この時、後で詳細は説明するが、シリコンと結ぶ直線が
[110]Siに平行になる為には、ペロブスカイトの
C軸が傾く必要がある。こうすることによって、本発明
者らの計算によれば、単に,ペロブスカイトとSiを結
ぶ線を並行に配置させて、C軸を傾けない場合にくら
べ、全エネルギがさらに10%低下することをも見い出
した。この上で更に、実際素子への利用を考え、PZT
中に0.01%以内までで有れば、充分単結晶ペロブス
カイトの本来の特徴を示す事を見い出した。
【0290】ここで指摘する界面エネルギの問題は以下
に説明するように、素子特性には多々影響を及ぼすもの
である。強誘電体膜と基板とが形成する界面に発生する
不整合によるエネルギは、不対原子を形成しやすくし、
これにより、界面準位の形成を助長することになる。こ
の界面準位は、伝導現象の妨げや、信頼性の低下につな
がる。また膜中においても、異常ホ゛ント゛等の欠陥がで
き、準位になる可能性をはらんでいる。ここでは、これ
らの欠陥の許容範囲を初めて示すものである。
【0291】本発明によって作成した、MIS特性と従
来例との比較をも行った。ここでは、特に界面準位を効
果的に評価できるものとして、CV法を用いた。本発明
による特性には、データのばらつきも正確に示した。本
発明によった場合、その界面準位に値は従来例に比べ、
約1/10に低下していることがわかる。また本発明に
よるMIS素子の移動度も同一電界状況下で評価すると
約12%の向上が見られた。また信頼性も向上した。
【0292】本発明者等は今までにない新しい軸関係を
発明するに至った手順をさらに詳しく説明する。図35
が本発明者等が作成した、全く新しい厳密abinit
io/分子動力学シミュレータの構造図である。図36
の上半分は、abinitioの部分で、下半分は分子
動力学シミュレータの部分である。このシミュレータを
用いて、初めて、発明者らは、例えば ペロブスカイト
の構造が如何なるものかを示した。図37が、本発明者
らが、初めて厳密な構造を示すことができた、ペロブス
カイトの透視図である。図中の大きい玉は酸素を示し、
小さい玉は、Pbを示している。ここでZ軸はペロブス
カイトのC軸を示しており、図でわかる様に、ペロブス
カイトのどの面とSiのどの面が巧く接合するかは、こ
の様に複雑な面であるので、極めて推測は困難であるこ
とが良く分かる。
【0293】また、上記の分子動力学シミュレータの精
度チェックとして、酸化膜に着目し、分光データの対比
作業を進めた。即ち、時々刻々のSi原子や酸素原子の
位置等から、O-Si-OやSi-O-Siの、ロッキンモードやス
トレッチングモード等を読みとり、これらの数値から固
有周波数を求めた。これと実測値とを比べたところ、そ
れぞれ、440cm−1と1100cm−1に対し、490cm−1及び
1111cm−1であり、妥当な範囲である。また圧力を印加
した場合、O-Si-Oモードの強度が増大することも実測値
と良く一致していることが分かった。
【0294】これらの検証作業をもすすめながら、本発
明者らは、PZTに非常に大きな圧縮応力を印加した場
合についても計算してみた。界面としては、この様な構
造をとったほうが、全エネルギの利得があることもわか
った。
【0295】計算に使った領域等について若干詳しく論
じてみる。計算ではSi基板側に、深さ方向に数十原子
層をとり、また、表面での奥行き方向に数十原子層をと
り、また、左右方向にも数十原子層をとり、計算領域を
とった。またペロブスカイトにおいても、縦横高さをそ
れぞれ数十原子層ずつとった。
【0296】今までにない手法で、新しい、位置関係を
見いだせたのは、この方法で摘めている。本発明者ら
は、強誘電体膜及びSi単結晶を含む領域に対して、厳
密なイオン結合ポテンシャルの厳密積算式を開発した。
従来の手法では、厳密性に欠け、特に系のエネルギを正
確に求めることが出来なかった。本発明者はこの明細書
のなかで従来の欠点や不正確な場所に付いて順次明らか
にして行く。
【0297】本発明者によると、以下のように説明でき
る。本発明者らは、従来にない新しいシミュレータを構
築し、一例として、強誘電体膜の単結晶、具体的にはペ
ロブスカイトについて計算を行った。
【0298】以下に、本発明者が作成した、シミュレー
タについても従来との比較をしながら、若干説明を加え
る。たとえば、まず、酸化膜について述べてみる。Si
−O間、O−O間、Si−Si間の3種類のポテンシャ
ルを厳密に表現しなければならない。Si−O間とO−
O間、Si−Si間のtotalのポテンシャルは実際はす
こぶる難しく、3種とも距離rの項により積算量は、無
限大に発散することになる。また計算自体もクーロンポ
テンシャルは遠くまで裾を引くので打ち切ることができ
ず、少し厳密な計算では従来の手法ではEwaldの方法を
使っていた。しかし、後述するが、計算に厳密性を欠い
ている。
【0299】実際に配置するに当たり本発明者らが用い
た手法を以下に示す。ここでは、本発明者らが初めて作
成したシミュレータを以下に示す。
【0300】本発明者らは、鋭意検討し、計算物理学に
則りながらも、新しい独自のabinitio分子動力
学システムを構築した。この新しいシステムがあったか
らこそ、新しい試みが可能になった訳である。以下に本
発明者等が作成したものを詳解する。
【0301】ここではSiO2を例に採っているが、他
のものについては、単に元素の番号等を入れ換えるだけ
で使用できることも確認した。
【0302】本発明者等が提案した新しいシミュレータ
の全体構成を図38に記す。分子動力学部分(MD)と
密度汎関数(DFT)とを合成した新しい手法を提案し
た。まず本発明者らによる分子動力学法と密度汎関数法
との使い分けを図39を用いながら以下に示す。
【0303】図36に示す様に、分子動力学(MD)部
分はフローチャートの下半分を占め、密度汎関数部分
(DFT)は上半分を占めている。分子動力学法は零K
でない所謂有限温度でのやや大きな系を扱い、原子(イ
オン)間のポテンシャルには、予め第一原理から求めた
ポテンシャルを用いる。他方、密度汎関数(DFT)
は、小さな系の基底状態を求める時に使うことにした。
そして、当然DFTの部分は温度は零Kとした。ポテン
シャル部分は、電子系の基底状態から求めた。具体的に
は波動関数(KS方程式)を解き、局所密度汎関数を併
用した。
【0304】本発明者等によるシミュレータの計算手続
きを実際にSiとOを例にとり説明してみる。DFT部
分で、縮重のない電子系の基底状態の全エネルギを求め
て置
【外1】 ]である。
【0305】即ち
【数104】
【外2】 ン相互作用エネルギー、第3項は、電子の1体近似を想
定して、電子間相互作用のない時の運動エネルギーであ
り、
【数105】 と書かれる。第4項は交換相関エネルギーで、その形は
局所密度近似による。電
【外3】
【数106】
【外4】 方程式(Kohn-Sham equation)
【数107】
【外5】 を自己無撞着(self-consistant) に解かなければならな
い。
【0306】
【外6】 ルギーεxc(n)に電子密度nを掛けて積分して、
【数108】 とした。εxc(n)の形はWigner他の式など、いろいろ提
案されている。さらに、
【外7】 のクーロン場をそのまま使うのでなく、核の位置での特
異性を消したノルム保存擬ポテンシャルでおきかえる。
このようにならすことによって、フーリエ展開での成分
の個数をおさえることができた。
【0307】全エネルギーとして、本発明者等は厳密を
期するため、原子(イオン)間のクーロン・ポテンシャ
ルも含むものを使う。即ち、
【数109】
【外8】 ε など)。式(2071)から右辺第1項はDFT
の式(2065)から右
【外9】 である。その部分はDFTの場合と同様に1体近似の式
(2065)からVeff
【外10】 ここで新しい考え方を入れた。Eをポテンシャル・エネ
ルギーと見なして、他に仮想的な運動エネルギー
【数110】 を導入する。従ってラグランジアンは L=K−E …(2073) 式(2073)中のMIは本当の原子(イオン)の質量
だがμ、μ は仮想的な質量であり、後述のように目的
に応じて変えた。束縛条件として正規直交条件
【数111】 を課するから、オイラー方程式を作る時に未定乗数Λik
を導入して、Lに
【数112】 を加えたものの変分をとる。その結果
【数113】 式(2073)の運動エネルギーの温度Tが対応する。
古典統計力学のエネルギー等分配則によって各自由度が
温度を持っているのだから、Σ1/2MII 2との関係でい
えば仮想的でなく現実の温度である。またμ、μ の大
小によってνiνを抑制したり増強したりできる。
【0308】例えば、μ MIとして高温からT→0と
すれば、RIをとめたままψiが変化して電子系の基底状
態を得る。この時式(2076)の左辺が0になり、右
辺は式(2071)とその下の注意により
【数114】 となるからψiの線形結合を適当に作る{Λik}が対角
線化されてKS方程式(2064)を得る。つまりDF
TでKS方程式を解いたのと同じ結果をsimulated anne
aling で得たことになる。ただし、温度Tの状態を作る
のに、モンテカルロ法によらず、仮想的な力学系の運動
を追っているから、dynamical simulatedannealing(DS
A)と呼ばれる。
【0309】式(2076)の積分をする時、必ずしも
ψiをそのまま変数として扱うのでなく、それを平面波
で展開した展開係数(つまりフーリエ係数)を変数とし
て扱うことができる。この展開係数の個数をむやみにふ
やさないためには、DFTの場合と同様に擬ポテンシャ
ルで核の位置での特異性を消しておくことが必要であ
る。特に結晶の周期性を利用できるときには、上記の
{Λik}を対角線化するψiの線形結合をφnkと記す。
kはブリユアン帯域内の波数ベクトル、nは帯域インデ
クス。
【0310】φnk
【数115】 を満たすはずである。λnkは行列{Λik}の固有値であ
り、エネルギー準位である。φnk(n、t)を展開して
【数116】
【外11】 満たされない)有効ポテンシャルVeffも展開して
【数117】 式(2075)と式(2076)を式(2074)に代
入して
【数118】
【外12】 共通になるから外して
【数119】 を得る。この式は振動部分の積分を解析的にすましてい
るから、純粋に数値的な方法よりもΔtを大きくして、
計算量を減少することが出来る。
【0311】クーロンポテンシャルは、以下の様にな
る。本発明者は分解してゆくと4項に分かれる事を見い
出した。特に従来は3項しかなかったが、新しく第4項
があることを確認した。これらの方法を順次しるす。こ
れらは以下に示す様に、本発明者らは、まず厳密な式の
出発として、誘電率を加味した基本式から解き始めた。
まず基本式である。
【0312】
【数120】 導体ε=∞と真空ε=1の場合とでクーロンポテンシャ
ルが異なり、Lは(立方体の)単位結晶の一稜、Σは単
位結晶内でとり、Zi、酸素i第i粒子の荷電と位置で
ある。これは球内の荷電によって球の内面に分極が生じ
ることによる。導体でない球の内面に双極子の層が出来
るが、上記式の最終項がちょうどその効果を打ち消す働
きをする。Ewaldの方法は左辺、ε=∞のものを与える
から、真空内の値を求めるには上記式の最終項を加えな
ければいけない。結果だけ掲げる。この部分でも従来の
導出はしばしば誤りがある。
【0313】クーロン・ポテンシャルを
【数121】
【外13】 、以下のように設定した。
【0314】
【数122】 と表される。ここに、nx,ny,nzがそれぞれx,
y,z方向の稜のベクトルで、nx,ny,nzはそれぞ
れ(バルク結晶の場合)−∞から+∞にわたる整数であ
る。Σ′の´はn=oの時のj=iを除くことを意味するも
のとする。
【0315】ここで新しいF関数を導入して、
【数123】
【外14】 を見い出した。
【0316】これを本発明者等はさらに下の様に変形し
て、
【数124】 指数mx,my,mzはそれぞれ−∞から∞にわたる整数
である。 = の時は、
【外15】
【数125】
【外16】
【数126】 を使った。
【0317】これにより、始めのクーロンポテンシャル
は、
【数127】
【外17】
【数128】 これは を含まないから力には関係ない。さらに次ぎの
ように変形できる。
【0318】
【数129】 今までの結果をまとめて、クーロン・ポテンシャルは、
結局、
【数130】
【外18】 /Lx,0,0)+my(0,1/Ly,0)+mz(0,0,1/Lz)である。上記の
表現の具合のいい点は、もとのΣの項が逆数のオーダー
でしか減衰しないのに対して、Φ(1)ではΣの項がer
fcの因子によって、Φ(2)ではΣの項がexpの因子に
よって、急速に減衰することである。Φ(3)での減衰
との遅速に逆向きに効くから、両者のバランスがとれる
ような適当なκを選ぶ。距離の近いところから順にクー
ロン力の寄与を計算した結果であり、しかも周囲が導体
である場合の結果である。周囲が真空である場合によれ
ばもう1項が加わる。これが特に従来省みられなかった
部分である。
【0319】本発明者等が作成した回転点群を図40を
用いて以下に示す。Si原子を半径rの球の中心(0,0,0)
に置き、球面上正四面体の頂点の位置に4個の酸素原子
を置く。C原子の位置は、オイラー角(θ、φ)で指定で
きる。それを(x,y,z)座標で表すには、北極に向かうベ
クトル(0、0、γ)に対して、まずy軸まわりθの回転、つ
いでz軸まわりφの回転をしてやればよい。
【0320】その回転行列は、
【数131】 ここで、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)に作用させて、
【数132】 となる。北極を指していた第1の酸素原子を(θ、φ)の
向きにするには、式(2061)のR(φ、θ)をそのま
ま作用させればよい。その結果、本発明者等は以下のよ
うに求めた。4個の酸素原子の(x,y,z)座標は、
【数133】 となる。単位格子内の4個のSi原子の(x,y,z)座標をパ
ラメタa,b,cで表す。図40と照合すると、2aが縦方
向、横方向に共通の周期である。
【0321】 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)の平行移動 を作用させれば得られる。
【0322】その他に、I II1とIII2、II III1とIV2、I
II IV1とV2の組がそれぞれ同一の酸素原子であるが、こ
れらは(当然のことだが)上と同じ関係式を与える。酸
素原子I1とII4は、横方向に1周期分2aだけずれてい
る。
【0323】この様に調べてきた上で、さらに実際プロ
セスを考えてどの程度の欠陥密度まで許容されるかにつ
いて調べた。即ち、今、例えば3200個のSiと32
00個のTiを用意し、これに0個から100個の酸素
欠陥を作った。この場合の時々刻々のSi原子や酸素原
子の位置等から、次の事が明らかに成って来た。まず、
酸素原子及びTi原子の動的挙動とPZT構造因子との
関連を調べてみた。単結晶PZT中に酸素欠陥が存在す
ると、数原子先の遠方のTi及び酸素までも、無欠陥の
場合と比較すると明らかに擾乱を受けており、しかも、
点欠陥の移動度はすこぶる大きいことである。これらの
動的挙動を、単結晶Siの場合と比べると、PZTの場
合の方が極めて影響も範囲も大であることが分かった。
この様な活発な動きは、PZTにはイオン結合成分がか
なり有り、分極しやすい性格を持っており、従って、局
所的な欠損が引き金に成ってポテンシャル変化に連動す
るものとも考えられる。本発明者らは、このような手続
きで、詳細に調べたところ、0.01%以下であれば充
分単結晶強誘電体膜の利点を引き出すことが出来ること
を見い出した。
【0324】本発明者による強誘電体材料の原子レベル
計算について以下にします。本発明者等は、PTOやP
ZO、PZTの原子レベルでの構造把握や、構造と膜特
性との把握計算を進めた。計算の結果を報告してみる。
原子レベル計算で用いている対象構造としては、(1)P
TOの完全結晶で、酸素、Tiがそれぞれ高対称位置に
あるもの、(2)PTOの完全結晶で、酸素が、0.3Aずれ
たもの、(3)前記(2)に酸素欠損を、c軸にそった部分に
1個入れたもの(図37)、(4)前記(2)に酸素欠損を、
a軸にそった部分に1個入れたもの(図43)、(5)P
TOのTをZrに置き換えたもの等である。特に、(5)
においては、Zrを、人工的に層状に置換したものも実
行した。酸素やPbの働きは、PTOやPZO等で異な
っていることが初めて分かった。これらから、強誘電体
膜の「自発分極」のきっかけが掴めた。また、強誘電体
膜の原子配置や歪み、欠損構造等と、基本膜特性との関
連もつかめた。
【0325】最も基本的なPbTiO3膜の原子レベル
での計算についてまず示す。引っ張り歪みを受けた場合
や、試みにPb位置を平衡点からずらせた場合(図4
4、図45)や、酸素欠損を有した場合等について、膜
特性を求める上で最も基本になる電荷分布が上手く求ま
るかを調べた。図44はPbTiO3の単位胞の絵で、
さらにPbを中心位置から、少しずらせた場合である。
この時の電荷分布の等濃度図を求めたものである。また
図45はやはり単位胞であるが、Ti-酸素結合の部分を
ひずませた場合である。この結果から以下の点が分かっ
てきた。(i)全体に非常にイオン結合成分が多いこと。
また、(ii)電気陰性度の強い酸素に殆ど大部分の電荷が
引き寄せられ、次にTiの廻りに電荷があることが分かっ
た。また(iii)このままでは電荷分布が対称的で自発分
極は起こり得ないが、酸素「欠損」や、Tiの位置変動を
考えると、局所的電気的中性が破られることが、容易に
わかる。さらに引っ張りにより、誘電率が低下すること
も確認出来た。
【0326】さらに例えば、25℃で熱擾乱がないとし
て見たときの計算結果では、以下のことが分かった。即
ち、点群C4vからなる理想的なPbTiO3膜の安定構
造を求めた。その結果、酸素を極わずかずらしても、エ
ネルギの安定な元の点群C4vに戻ることが分かった。さ
らに大きく、酸素位置を0.3Åずらしたところにも、通
常言われている様に、もっと安定な平衡点があるか。こ
れを調べた。出発点0.3Åより極わずか離れた位置に向
かってさらに安定になりつつある。この辺りが通常言わ
れている位置で、このずれ故に自発分極が生じるとも言
われている。安定点が少なくとも二つあり、1つは単純
立方晶の位置にあるもの、もう1つは若干ずれた位置で
あることが初めて分かった。
【0327】分極の起源は、局所的に電荷量の「中性」
を、くずす様な原子位置が、自由エネルギ的に安定点で
あることとも置き換えられる。また、以上の計算では、
酸素の存在は極めて大きいことが分かった。それ故、酸
素欠損は、ひょっとしたら、結晶的にも、電子論的にも
大きなトリガーになりうることがわかる。即ち、 (1)各3種類の原子にそれぞれ電荷分布は局在している
様である。
【0328】(2)周期律表にある順序に従って、O及び
Ti及びPbは、それぞれ(1s2,2s2,2p4),(1s2,2s2,2
p6,3s2,3p6,3d2),(1s2,2s2,2p6,3s2,3p6,3d10,4s2,4p6,
4d10,4f14,5s2,5p6,5d10,6s2,6p2)である。これを反映
して、特にPbの電荷分布領域が大きいのが分かる。
【0329】(3)良く見ると、とりわけ、Tiの電荷分
布は直ぐ隣のOによって変調されているのが見える。
【0330】これを抑止するため、F,Cl、Br、I
等のハロゲン元素を上手く置き換えてやる手も有ること
が分かった。
【0331】例2 上の例では、Si・PZTについて
記したが、単結晶強誘電体膜として、スピネルMgAl
24膜、酸化セリウムCeO2 、チタン酸ストロンチウ
ムSrTiO3 を用いることもできる。また、酸化アル
ミニウムAl23、酸化ジルコニウムZrO2 、酸化イ
ットリウムY23、イットリウム安定化ジルコニウムY
SZ、PrO2 、弗化カルシウムCaF2 など、およ
び、その積層膜をも用いることができる。また、該単結
晶強誘電体膜と下地シリコンウエハあるいは下地電極と
の界面にシリコンPZTを挟む変形は本例以外の構造に
おいても応用することは可能である。
【0332】また、電圧下での自由エネルギを小さくす
ることが大切であり、歪のエネルギだけに着目すれば、
そのときは必ずしも最低になっていなくても良い。我々
がペロブスカイトについて示したのも、思想に基づいて
いる。
【0333】本発明による実際の測定系についても説明
する。図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できまる。損失のない常誘電体
をこの回路で観測すると、原点を通る直線になるが、損
失のある場合や、強誘電体はヒステレシスをもつ。図4
8に本発明になるD−E曲線を実線で表し、比較のため
に従来のPZTの結果を破線で示した。これから分かる
様に、本発明になるものは、非常に誘電率が大きいこと
もわかる。
【0334】ここで指摘するところは、強誘電体膜と基
板とが形成する界面に発生する不整合によるエネルギを
印加時で最低にしようとするものである。しかも、この
条件で実際利用上のことを考えその品質の許容範囲も示
した。この界面エネルギを縮小することは、伝導現象の
妨げになる準位の形成抑止や、信頼性の向上につなが
る。この構造を本発明によれば、新しいabiniti
o/分子動力学シミュレータによって予測するものであ
る。所望の評価関数を搭載したabinitio/分子
動力学シミュレータは、経験的なモデル等を一切用いて
いないので、このシステムに指摘するところは、全て自
然現象を完全に予測できるものである。ここでは、系の
印加時の全エネルギに着目した例を示しているが、例え
ば誘電率を所望の値にすべく設計についても、さらに
は、劣化抑止を所望の問題であれば、それに適した設計
が可能になる。
【0335】以上2次元の物性を計算する方法を示し
た。3次元の物性の計算方法の一例については、図8と
共に既に説明した。
【0336】叉、実際の対比や診断がどの様にして行わ
れているかをしめそう。即ち、シミュレーション結果
は、スケジューラ(13)に転送される。この時、他
方、実際の実験系からは、モニターデータが、時々刻々
送られ来る。そして、これらは、ファイル(4)に格納
されるととものに、ファイル(4、3)からのデータ転
送は、スケジューラによって同期させておく。そして、
対比診断検定部IVにある対比・検定・診断ブロックに
転送される。
【0337】図49は、本発明に使用されるクラスタ化
ツールにおいて、in−situモニタ測定装置がどの
ようにして付加されているかを示す概念図である。図4
9に示すように、モニタ量は、例えば、入力ガス構成分
析やXPSシグナル、FT−IRシグナル等であり、こ
れらは、所定の時間空間内における平均値であり、ほん
の一部分の情報である。本発明者らは、これら限られた
測定量と、他方計算から得られた十分な情報とを対比さ
せる方法を提供するものである。これらは、図50の様
に、モニター上のグラフィックな画像として表示され
る。
【0338】図51は、本発明によるabinitio
分子動力学シミュレーションシステムへの入力データの
一例としての温度シーケンスを示すものである。また図
52は、本発明によるabinitio分子動力学シミ
ュレーションシステムによる、原子レベル実行予測の概
念を示すものである。このように、例えば原子レベル
で、基板上での反応や、分子の解離などが記されてい
る。
【0339】図53は、本発明によるabinitio
分子動力学シミュレーションシステムにより得られた原
子レベル実行予測の可視化図である。これは、非晶質S
iが堆積されたあと、600℃近傍の高温でアニールさ
れた状態を原子レベルでみたものである。
【0340】また、本発明のabinitio分子動力
学シミュレーションシステムは、Siのみならず、酸化
膜に対しても適用することができる。即ち、図54は、
本発明によるabinitio分子動力学シミュレーシ
ョンシステムのSiO2についての実行例を示す。ここ
では、SiO2におけるSi−O−Si角が時々刻々求
められており、これから、フーリエ変換し、固有振動数
をもとめ、これを重畳することにより、FT−IRのス
ペクトルが計算のみで求められる。この予測結果を用い
て、実測値の時々刻々データと比較検定すれば良い。
【0341】図55は、さきに示した図51の入力デー
タを用いた場合の、本発明のシミュレーションシステム
を用いた原子レベルでの実行例を示す。図55から、予
想どおり実行されていることがわかる。
【0342】図56から図69は、図51の入力データ
を用いた場合のSi膜の成膜から、アニール、さらに固
相成長までにおけるSi−Si間距離、結合角を、時間
と温度を変化させてもとめたものである。
【0343】図1に示したシステムでは、基本膜特性の
算出を膜特性算出ユニットIIIで行った。しかし、図
70に示したシステムのように、特別な膜特性算出ユニ
ットIIIを省略し、プロセス変動算出ユニットI内部
で必要な特性の算出を行ってもよい。それらは設計の要
請に応じて、適宜変形可能である。
【0344】以上のように、本発明によって、ハードウ
エアや計算量を従来に比較して大幅に軽減することが可
能となった。例えば、計算量について言えば、プロセス
シミュレータTOPAZでは、図98に示すように、メ
モリは約24Mバイトで、実行時間は数分間となる。こ
れに対して、有名なプロセスシミュレータSUPREM
4では、粘弾性等を考慮しているので、メモリは200
Mバイト以上、実行時間はクレイで1日にも達する。さ
らに分子動力学やab−initio計算になると、メ
モリが数倍、実行時間も数倍にも達している。
【0345】このような状況では、分子動力学やab−
initioツールをフルに用いて、原子の動的挙動の
予測や、光学的諸特性の予測には、計算機のコストが大
幅にかかかってしまう。本発明者らは、これら不利益
を、新しい数学的な技巧を駆使し、克服する方法を見出
すとともに、この技巧により、シミュレータの新たな利
用技術をも見出したものである。
【0346】従来でも、上記に若干の手続きを加えたも
のもある。それが図99に示すフローシートであり、従
来の分子動力学シミュレータにおける密度汎関数の流れ
の主要部分を示す。即ち、この手法では、仮想質量を与
え、全系のラグランジアンを用いる。これは、本発明者
等が先に記した式(1223)に該当する。これはすで
に知られているCar−Parrinello法に相当
するものである。この方法は一見高速に計算できるとさ
れている。しかし、それでも例えば、クレイ社レベルの
スーパーコンピュータでも、1つの場合の計算でも1昼
夜を要する。従って、種々のパラメータの変動を調べる
ため、1回1回計算をしていたのでは、到底、実時間で
の診断ツールとはなり得ない。
【0347】又、プロセス装置を実際に稼働させるに
は、それなりの準備と費用がかかる。本発明によれば、
半導体製造システムを構築し、実際に半導体の作成プロ
セスを進行させることなく、プロセスシミュレーション
あるいは思考実験のみを行なうことができ、これは非常
に有用である。例えば、様々な分子構造モデルを用いて
思考実験を行ない、それを用いたデバイスの特性を予測
することによって、新たな材料の探索を手軽でしかも効
果的に行なうことを可能となる。
【0348】図131はそのようなプロセスシミュレー
タの構成図であり、コンピュータ1、ハードディスク装
置などの内部記憶手段3、マウス5やキーボード6など
の入力手段、モニタ7やプリンタ9などの出力手段、及
びCDーROM11やフロッピーディスク13等の外部
記憶媒体及びその駆動ドライブ15からなっている。
【0349】内部記憶手段3には、予め分子動力学シミ
ュレーションを行なう為の手続きを記述したプログラム
が作成され格納されている。このプログラムは、例え
ば、フォートランやC言語等の高級言語に必要によりア
センブラ言語を組み合わせて作成される。勿論、CDー
ROMやフロッピーディスク等の外部記憶媒体に予め保
存されているそのようなプログラムを、適宜内部記憶手
段3にインストールしてもよい。
【0350】ここで行われる処理は、すでに説明したも
のと同じであるが、実際に稼働している半導体製造シス
テムとのインターフェイスは省略される。つまり、図1
の分子動力学シミュレータII、ファイル1、原子レベ
ルまで考慮しないマクロ的な性質を計算を行なうプロセ
ス変動算出ユニットI及びファイル2を含むループに対
応する処理のみが行われる。
【0351】又、この様なシステムは、図1の実際的な
半導体製造システムのプログラムから対応部分を切り出
して利用することが可能である。又、逆に、このような
計算のみのプロセスシミュレータをコンピュータ上で実
現、実用化した後、実際のプロセス装置との組み合わせ
を行なうことも有り得る。即ち、図1の半導体製造シス
テムは、プロセス装置V、スペクトルメータVI及び図
131のシステムから構築することができ、スペクトル
メータVIの外部出力端子を、適当なケーブル19を用
いてコンピュータ1の通信ポートへ接続して、実測値と
計算値との比較処理が行われる。その結果は、同様に適
当なケーブル21を用いてコンピュータ1の通信ポート
へ接続されているプロセス装置の制御手段23へ送ら
れ、リアルタイムでプロセス制御に反映される。
【0352】又、酸化やシリサイド等の工程誘起応力に
関しては、すでに2次元シミュレータは、ほぼ実用化レ
ベルまでに達している。しかし、以下に説明するよう
に、本件発明者は微細素子の設計には3次元が必須であ
るということを発見し、本発明を3次元シミュレーショ
ンとして完成させた。
【0353】今までのシミュレーションでは、酸化にし
ても拡散にしても、結局のところは、Fickの拡散法則に
帰属する。x‐y面を扱う2次元問題では、第3軸である
z方向には物質の移動がないとしてよいが、「応力」の
場合は、物質が動くのではなく、状況が大きく異なるこ
とが分かつた。これは、酸化に眼らず、シリサイデーシ
ョンでも同じである。
【0354】即ち、現状でもつとも進んでいると言われ
ているTSUPEM4の2次元での応力計算は、応力を
σとし、変位をvとし、粘性をμ、ボアソン比をνとし
て、以下の式を基本としていた。
【0355】
【数134】 σxx+σyy=μ/(0.5−ν)・(∂vxノ∂x+∂vy/∂y) σxx−σyy= 2μ・(∂vxノ∂x−∂vy/∂y) σxy= 2μ・(∂vyノ∂x+∂vx/∂y) これで分かることは、2次元的な単純なxy面の取り扱
いで、ここではz成分が出現していないことである。ま
た、ここでは陽に見えていないがヤング率Eはスカラー
量として扱っている。
【0356】本件発明者は、厳密3次元剛性マトリック
スの作成と厳密面方位演算子の予備計算を行なった。そ
の結果だけを記すとSi系では,以下の様になることが
分った。すなわち、
【数135】 この式を用いて、まず、奥行き方向に十分厚い2次元的
平面歪み(plane strain)の意昧を考えてみる。即ち、
(εx≠O、εy≠O、εz=O、εyz=O、εzx=O、εxy≠
O)の場合について、上記式に代入するとσzの項が残存
るのが分かる。すなわち、σz=Cxz・εx + Cxz・εy
+ Czi・εxy。またσzの値は、決して微少量ではない
ことも突き止めた。従って、3次元は必須であり、上記
TSUPEM4では、2次元でも考慮すべきσzの項が
欠落していることが初めて分かつた。2次元であつて
も、少なくとも、酸化やシリサイドの応力の場合は、第
3軸のz成分も出てくることが明らかになった。
【0357】次に、上記式を変形し、実際の各面での剛
性マトリックス Cijを求めてみた。単結晶Siの弾性マ
トリックスをまず記述する。
【0358】
【数136】 であり、C11=1 7 2±15Gpa、C22=66±6G
Pa、C44=−82±9GPaである。これで分かるこ
とは決して単一のスカラ―量ではないことである。これ
に、面回転演算子を操作させ、まず通常よく見かけるシ
リコン結晶の断面(110)について求めると、以下の
様になる(図132(B)参照)。
【0359】
【数137】 これで分かることは、z成分の寄与がきわめて大きい事
が分かる。また、バイポーラ等で問題になる(111)
面上については以下の様になる(図132(B)参
照)。
【0360】
【数138】 このような、変換操作は必須であるにもかかわらず、こ
れまでのシミュレーションでは考慮されておらず、微細
構造の解析の精度が低いということが判明した。そこ
で、本発明では、上記のごとく3次元を取り扱うことと
なった。これが、本発明の別の様相である。
【0361】
【発明の実施の形態】以下、本発明の種々の実施例につ
いて説明する。
【0362】実施例1 この実施例は、本発明の大きな特徴の1つである光学的
スペクトルの直接予測機能に係るものである。ここでも
う一度、本発明者らが考え出した手法から詳細に記述し
てみる。
【0363】まず、高信頼化酸化膜設計に関する技術お
よびその評価関連技術等について、従来の技術的レベル
から少し記述してみる。公知刊行物としては、例えば
「ThePhysics and Technology of Amorphous SiO
2 」Plenum Publishing Corporation,ISBN 0-306-42929
-2がある。この刊行物には、種々の酸化膜の光学的測定
結果として、IRの結果やラマンの結果が記載されてい
る。この刊行物は、特に、原子レベルの議論としては優
れたものと言われている。
【0364】この刊行物のp.63には、ポール・マク
ミリアン(Paul McMillan )によって非晶質酸化膜のラ
マン信号が、工程によってどの様に変化するかが論じら
れている。図71及び図72は、ポール・マクミリアン
によって測定された結果である。例えば、図71の曲線
aは、通常の酸化膜のラマン測定結果を示し、曲線b
は、デンシファイした酸化膜のラマン測定結果を示す。
【0365】従来では、この両者のスペクトルを比較
し、例えば、デンシファイにより水分が脱離し、これに
よってスペクトルが変化したものとされていた。また図
72は、単結晶水晶に圧力印加する前後でのスペクトル
をみたものである。特に610と500cm-1のところ
に変化があるとしている。しかし、これでわかるよう
に、まず一見して、スペクトルの変化量はきわめて小さ
く、どこに着目すべきか、またそれを有為差とみるべき
かについては大変技術を要するものである。またさら
に、そのスペクトル変化量から物質の変化量を同定する
にはさらに複雑な経験を要するものである。
【0366】この様な状況なので、今後さらに必要にな
ってくる新しい材料の開発においては、スペクトルの同
定テーブルがないものが多々あり、スペクトルの解釈が
大きな障害になってくる。
【0367】また、上記刊行物のp.71からは、小林
等が高応力下での酸化物ファイバーのラマンスペクトル
の変化をみている。その結果の一部を図73に示す。こ
こではSiO2ファイバに応力を無印加の場合と3.5
%の引っ張りひずみを印加した場合について示してい
る。図73から分かるように、3.5%の引っ張り歪み
を加えた場合、450cm-1近傍のところに変化があ
る。しかし、ここでも分かる様に、スペクトルの変化量
はきわめて小さく、やはりどこに着目すべきか、またそ
れを有為差とみるべきかについては、大変技術を要す
る。また、そのスペクトルの変化量から物質の変化量を
同定するにはさらに経験を要する。
【0368】上記刊行物において小林は、図74を用い
てそのスペクトルの解釈を試みている。図74から分か
る様に、引っ張りひずみを受けるとSiO2構造の10
9,5度の4面体部分が歪むとしている。しかし、彼
は、たった1組の4面体構造を持ち出しているにすぎ
ず、となりの4面体部分とはどのようなつながりになっ
ているか全く議論していない。実際には、ある有限の温
度であるので、必ず角度の分布を伴うはずである。
【0369】本発明者がこの明細書で何度も記述してき
たように、定量的な分布を把握しないと、およそスペク
トルの解析には程遠いものとなる。また図33も同じ議
論が小林によって展開されている。図33(a)は、O
−Si−O結合角の関数としてのSiO4 四面体の通常
モード振動の波長を示す特性図である。また、図33
(b)、(c)は、ストレスにより生じたV4モードの
デカップリングによる、強度の減少を示す特性図であ
る。
【0370】図76は、上記刊行物のp.107に記載
されたW.B.Fowlerによる酸化膜中の欠陥モデル図であ
る。図76の(a)は、完全α型クォーツの一部の結合
モデルを示し、図中Lは長結合、Sは短結合である。
(b)はE1 ´中心のモデル、(c)はE4 ´モデル、
(d)は、E2 ´モデルをそれぞれ示す。(b)〜
(d)における不対電子は、e- で表されている。
【0371】図76からも分かるように、上記刊行物で
は、109.5度の4面体構造を2個または1個用意し
て欠損構造を議論している。しかし、上述したように、
実際には、ある有限の温度で必ず角度の分布を伴うはず
である。従って、ある種の多数個の計算を行わなけれ
ば、およそスペクトルの解析には程遠い。
【0372】図77は、本発明者等が作成した新しい分
子動力学シミュレータによる実行例である。Siは32
個、酸素は64個を用いた。また逆格子点は520点で
あった。縦軸は体積の変化率、横軸は時間をそれぞれ示
す。図77から、1000ft後でほぼ外圧と平衡にな
ることが確められた。例えば、このときのO−Si−O
角、Si−O−Si角、及びO−Si距離を図24に示
す。
【0373】図24から分かる様に、各原子の「戸籍
簿」に相当する膨大データが必要なのである。本発明者
等は、計算手法にも幾つかの工夫を用意している。例え
ば、赤外吸収では、時々刻々の原子の振動運動の時系列
データから、赤外吸収の選択律を満たす振動モードを抽
出し、その頻度を求めればよい。いずれにしてもこの様
な試みはまだ発表がなく、これはプロセス開発の新しい
「兵器」と言うことが出来る。
【0374】本発明者等は、分子動力学システムを用い
て、十分、平衡に達したところで、時刻を0に設定し、
その後、時刻tまでを今サンプリング時間とする。時間
相関は以下の通り求めた。即ち、 である。ここでM(t)は分極率であり、下記式で表さ
れる。なお、上記式の積分値は−∞から∞まで積分した
値である。
【0375】 M(t)=Σei・r(t) …(1362) また<M(t)・M(0)>は、Mの自己相関関数を表
わす。
【0376】ここで、ei はその点における電荷数であ
り、iは粒子の番号である。またrはベクトルであり、
その測定起点は特に制約はないが、計算セルの角にとる
のがよい。
【0377】時間相関を求めるには、それぞれの時刻に
おけるs(t1 )・s(t2 )・s(t3 )・s(t4
)・s(t5 )・s(t6 )・s(t7 )・s(t8
)を求めればよい。
【0378】例えば、サンプリング時間全体に渡っての
時々刻々の座標位置の計算値やその運動のエネルギ、さ
らには運動量等が、全粒子に渡って存在している。しか
も、さらに互いの粒子間の距離や、それらの角度位置関
係等も計算されている。
【0379】本発明者等は、幾つかの結晶性酸化膜や、
また非晶質酸化膜についても計算した。その結果の一部
を図78に示す。即ち、図78は、1つの例として単結
晶クリストバライトのSi−O距離の計算結果を示す。
また、図79は、欠損がある場合のSi−O−Si角の
計算結果を示す。図79から、欠損の回りでは運動が大
きく乱れていることが分かる。これらの欠損の濃度を種
々変えて同様の挙動を求めることも行った。
【0380】図80は本発明者等が作成した新しい分子
動力学の実行例であり、原子の動きをシミュレートした
ものである。図中、大きい玉は酸素を示し、小さい玉は
Siを示している。図80の(a)は、XYZ空間にお
ける原子を示し、(b)は、そのXY面への投影を示
す。ここではその軌跡まで示せなかったが、この単結晶
例では原子は主にXY軸の中線方向に原子が優先的に運
動していることも分かった。このことは上記式に示した M(t)=Σei・r(t) …(1371) の部分に大きな偏りを示し、結果として、単結晶特有の
スペクトルを示すことになる。
【0381】本発明者らは、測定温度を変えて、個々原
子の挙動や、これを基にした固有振動の様子、さらに
は、歪みを加えた場合についても、シミュレーションを
実行した。図81及び図82はそれらの例を示す。即
ち、図81は本発明による新しい分子動力学シミュレー
タによるX線スペクトル予測図、図82は、同様の応力
予測結果を示す図である。
【0382】実施例2 本発明者らは、酸素欠損がある場合の酸化膜のスペクト
ルの予測も行えるようにし、実際の酸化膜との対比を可
能とした。実施例2として、酸素欠損酸化膜の光学的性
質まで計算可能とした例を示す。
【0383】ab−initio分子動力学法を使用し
た場合には問題にならないが、TTATポテンシャルや
BKSポテンシャルを使用した分子動力学法の場合に
は、酸素欠損を扱う場合に、各原子のチャージをどのよ
うに設定するのが最良であるかが問題となる。これは、
非常に歪みが大きくなると、結合が伸びたり、Si−S
i結合が生じ、一定電荷を仮定することができなくなる
からである。例えば、Si−O結合の伸びにより、実際
には電荷の分配が変化し、分極が大きくなると、伸縮振
動の静電エネルギーは大きくなり、IRのピーク振動数
は大きくなると予想される。
【0384】そこで本発明者らは、Rappe−God
dardにより、アルカリや水素を対象に試みられたこ
とのあるイオン化ポテンシャルと電気親和力を用いた平
衡電荷近似(Charge equiliburati
on approach)を、初めてSiO2に適用
し、電荷の補正を取り込んで静電ポテンシャルを計算
し、個々の原子の運動とIRを求める手法を行った。
【0385】なお、平衡電荷近似とは、系の静電エネル
ギーを各原子のイオン化エネルギーと各原子対の電気親
和力と各原子の電荷の関数と仮定して扱い、平衡状態で
は各原子の原子化学ポテンシャルが等しいとの条件か
ら、原子数に相当する数の連立方程式を解き、各原子の
電荷を求める手法である。
【0386】この方法をβ−石英に用いた場合のIRラ
インスペクトルと従来例を図83に示した。本方法を用
いることにより、ピーク周波数が約一割、低周波数側に
シフトした。この値は実測とよりよい精度で一致する。
a−SiO2でも本方法を用いることにより、より一
層、構造とIRスペクトルとが高精度で関連づけること
が可能となり、実測との対比が効果的に行えるようにな
った。
【0387】図84に、本発明者等が編み出した方法を
分かりやすく概念図で示した。また探索手法を図85に
示した。
【0388】なお、本発明は、上記実施例に限定され
ず。液相デバイスや超電導デバイス等、他のデバイスを
特にクラスタ化製造装置により製造することに対して有
効に適用可能である。その他、本発明の要旨を逸脱しな
い範囲で種々変形して実施可能である。
【0389】実施例3 ここでは、第3の実施例として、光学的スペクトルの直
接予測機能についての詳細を示す。はじめにa−SiO
2(アモルファスSiO2)の計算機上で作成し、その光
学スペクトルを予測する方法について示す。a−SiO
2は、単結晶SiO2を計算機上で溶融した後、急冷する
ことにより作成した。図86は、a−SiO2の作成の
際の温度設定を示す。溶融はゆっくりと行ってもよい
が、ここでは温度を急上昇させて急激に溶融させた。こ
の後、同程度の降温速度で降温した。図87はその際の
圧力の変化であり、図88は計算領域の体積変化であ
る。本発明者らは、この過程において、経験的分子動力
学法を用いた場合でも、ab−initio分子動力学
法を用いた場合でも、昇降温により温度を大きく変動さ
せる場合には、外圧と体積の扱いが重要であることを見
いだした。例えば、本実施例では、体積をほぼ一定に保
って相変態が起こるように外圧を温度の関数として扱
い、温度が上昇すると圧力が大きくなるように設定して
いる。その様子が図87,55に示されている。
【0390】しかし、体積に全く注意を払わず、外圧を
一定に保った場合には、図89に示したような体積の大
膨張が生じ、現実の高温溶融SiO2とは全くことな
り、著しく低密度のSiO2となる。更に、冷却過程で
も、固相に入っても有りうべからざる低密度のSiO2
が形成される。もちろん、ポテンシャルエネルギーを調
べると非常に高く、不安定な状態であった。こうなると
実測との対比は困難となってしまう。従って、外圧に十
分注意を払ってa−SiO2を作成することが望ましい
ことをここで見出し、実行した。
【0391】図90には、溶融時のポテンシャルエネル
ギーの変動を示した。外圧に注意を払っているので、ポ
テンシャルは振動しながらも緩やかに上昇している。急
冷の際にも同様に緩やかにポテンシャルエネルギーは減
少してゆく。図91にはこのときの計算領域のセルc軸
に垂直な面の2辺の時間変化を記した。セルは、相変態
に応じて形状が変えられるように、Parinello
−Rahmanの手法を採用しているので、溶融によっ
て変化している。
【0392】最初は石英が菱形であるため、2辺のなす
角度は60度となっている。ところが溶融すると、結晶
構造は壊れ、各辺に均等に内圧がかかるので、2辺のな
す角度は90度に近づいている。急冷しても石英には戻
らず、約90度で直方体を保つようであった。
【0393】更に、本発明者らは、a−SiO2を作成
する際、もう一点注意すべき点があることを見いだし
た。それは、SiとOのポテンシャルに関するものであ
る。ab−initioでも経験的ポテンシャル(AT
TAポテンシャルもしくはBKSポテンシャル)の分子
動力学法でも、SiO2系を扱う場合には、原子の配置
の変動によるポテンシャルエネルギーの変動が非常に大
きく、そのために、原子の運動を扱うには、通常のVe
lretのアルゴリズムよりも高精度に原子の運動を求
める必要があるということである。
【0394】図92は、通常のVerletのアルゴリ
ズムで原子の運動を求めた場合のポテンシャルエネルギ
ーの変動と運動エネルギーの変動、及び両エネルギーの
和を示す。外界とのエネルギーのやりとりが無い場合、
総エネルギーは一定になるべきである。しかし、通常の
Verletのアルゴリズムでは、ラグランジャンから
導かれる運動方程式の積分が十分な精度で行われないた
め、このようにエネルギーの保存が行われない。積分精
度を向上させるために、積分間隔である時間間隔を小さ
くとっても、状況はほぼ変わらなかった。これは調べて
みると、酸素のポテンシャルが非常に急峻であることに
起因する材料に特有のものであることが分かった。
【0395】そこで、酸素を正確に扱うために、Ver
letのアルゴリズムを改良し、時時刻々の個々原子の
運動の加速度の高次のまで求めて、これを積分して速度
や位置を求めて見た。すると、図93に示したように、
ポテンシャルエネルギーと逆位相で運動エネルギーが変
動し、総エネルギーは一定に保たれるようになった。こ
のようなエネルギーの正確な扱いがなされない場合、本
発明者らが実行した結果によると、高温で原子の運動が
激しい場合や、液体状のSiO2などで、計算精度の悪
さが原因で、構造が無秩序に壊れ、a−SiO2の作成
ができなかった。また、IRスペクトルも、全く実測と
異なるものになってしまった。
【0396】以上のように、本発明者らは、a−SiO
2の作成において独自に工夫をこらすことにより、初め
てa−SiO2を計算機上で作成することを可能とし
た。このようにして作成したa−SiO2は、石英に比
較して、Si−O結合長が若干長く、しかもSi−O長
は一定の範囲内に収まっていた。溶融SiO2では、S
i−O長は石英でのSi−O長よりも小さな結合長のS
i−Oが多く存在していたが、冷却するとこれが無くな
り、石英よりも長い結合ばかりになった。配位数を調べ
ると、第一近接を2Aとすると、Siは4配位が殆ど
で、3と5が少し混じり、それ以外は生じなかった。酸
素は2配位が殆どであった。
【0397】次に、a−SiO2に静水圧もしくは一軸
性の応力をかけ、種々の応力場に応じたa−SiO2
構造変化を予測計算するため、ポテンシャルエネルギー
が平衡に達するまで、温度を900℃に保って、2ナノ
秒の時間、SiO2を構造緩和させた。できあがった種
々のa−SiO2の構造を、個々、原子位置、速度、加
速度、その一次微分値、セル基本ベクトル、セル基本ベ
クトル変位速度、等の数値に変換し、それぞれ保存し
た。
【0398】ここまで行って、種々のa−SiO2構造
を作成した後、それぞれに対してIRスペクトルを求め
た。IRスペクトルは、ダイポールモーメントから求め
た。ダイポールモーメントMは、下記式 M=Σqi ・ri …(1461) から求めることができる。ここで、qiは原子iの電
荷、riは位置である。この位置を求めるとき、系が中
性の場合(Σqi =0)には、原点r0 を何処にとって
もダイポールモーメントの値は同じになる。それは、
【数139】 M=Σqi (ri −r0 )=Σqi ・ri −Σqi ・r0 =Σqi ・ri −r0 Σqi=Σqi ・ri …(1462) となるためである。また、系が非中性の場合には、r0
の掛かる項が残るので、原点の取り方によってダイポー
ルモーメントの絶対値は変化するが、フーリエ変換した
場合には周波数がゼロの成分に現れるだけで、他には影
響しないことが分かる。このダイポールモーメントを求
めるに当たり、原子位置rの取り方には一点注意が必要
なことに本発明者らは気づいた。それは、分子動力学法
を用いる際に周期境界条件を採用しているため、セルの
−x方向に出ていった原子が+x方向から入ってくるこ
とに起因するダイポールモーメントの不連続な変動を無
くすことである。
【0399】実際にIRを測定する温度では、a−Si
2は熱振動をしているのみであり、さほど大きく拡散
をしないことを考慮して、ダイポールモーメント算出の
際に用いるrは、周期境界条件を考えず、−x方向にセ
ルから出ていっても、セル外のその位置に電荷があると
扱って、モーメントを計算した。つまり、個々原子運動
計算の際には周期境界条件を考慮した原子位置と、周期
境界条件を考慮しない原子位置の2つの原子位置を両方
取り扱うことにした。このようにして求めたダイポール
モーメントの時間変化の例を図94に示した。3次元計
算であるので、x,y,z各方向のダイポールモーメン
トが得られた。
【0400】IRスペクトルを求めるには、第一段階と
して、このダイポールモーメントをフーリエ変換し、I
Rラインスペクトルであるパワースペクトルを求める。
その一例を図95に示した。この段階でIRスペクトル
のピーク位置が予測できる。すなわち、約1000cm
-1に大きなピークが出現すること、約500cm-1に大
きなピークが出現することである。この結果は、a−S
iO2のIR実測と一致した。
【0401】IRスペクトルを求めるには、このライン
スペクトルI(ω)にωと温度Tの関数を掛けて次の式
で求める。
【0402】
【数140】 α(ω)=(4π 2/3hcn)・ω・(1−exp(−hω/kT))・I(ω) …(1463) 種々の応力場でのa−SiO2のIRを求めることによ
り、応力場とa−SiO2の構造とIRスペクトルの詳
細(ピークシフト、ピークショウルダー等)が関連づけ
られた。これにより、実測との対比が可能になった。
【0403】ここで、図96に、本発明者らが新しく開
発した分子動力学シミュレータをα−SiO2に適用
し、計算されたIRスペクトルを示し、図65に本発明
者らが新しく開発した分子動力学シミュレータをa−S
iO2に適用し、外圧が加わった場合の予測されたIR
スペクトルを示す。また、図100及び図101に本発
明者らによる計算手順を示す。なお、図100における
ブロックB内の、個々の原子の運動の式におけるA、B
は、下記式に示される項である。
【0404】実施例4 本発明に係るシステムを用いることにより、工程中の誘
起応力の大きさやその成分の予測、さらには、その応力
と材料物性との詳細な関連の把握により、基板や層間膜
の塑性変形等の現象を予知することができる。この実施
例では、特に、応力問題を例にあげる。応力誘起をいか
に予測し、プロセス変動幅を考慮しながら、その低減を
はかるかを検討した。
【0405】昨今の微細素子プロセス設計には、プロセ
ス進行時の応力を考慮することが不可欠である。応力に
ついては、例えば以下の現象が生ずる。
【0406】(1)埋込み素子分離法では、Si基板に
溝を掘り、ここにSiO2を埋め込む。この場合、Si
基板とSiO2との熱膨張係数の違いにより、歪や応力
が残留し、熱履歴によっては、結晶欠陥を誘起する場合
もある。従って、できるだけ残留応力を低減できる素子
構造や製造プロセスを立案する必要がある。
【0407】(2)高品質の酸化膜を作成する目的か
ら、格子間酸素濃度の低いSi基板を用いることが頻繁
になってきた。このような状況では、Si基板強度が減
少し、塑性変形を誘起しやすい。
【0408】(3)浅い接合層の形成維持やウエル領域
の維持を目的に、アニールや酸化工程等の熱工程が、R
TPの採用により、一層基板に負担が多くなってきた。
【0409】このような状況のなかで、如何に欠陥の発
生を抑止するかとともに、工程途中で、モニタしなが
ら、応力が規定値やその変動予測値よりも大きい値が観
測されたなら、その原因を即座に明確にするとともに、
今後の工程において、どのように修正を加えるべきかを
示すことが重要である。このことは、クラスタ化ツール
では一層重要である。
【0410】さらにまた、応力の種類や同定は極めて重
要である。最大応力の主応力と主面を追跡することも大
切であるが、分解剪断応力に着目することも重要であ
る。例えば、Si(100)基板を用いている場合、そ
の最優先滑り面は図102に示すように(111)面で
ある。この(111)面上で[100]方向が滑り方向
であるが、この方向に剪断力が働いたときが最も危険で
あると考えられる。
【0411】本発明者らによるシミュレーションシステ
ムは、これらの同定の出力も可能である。このように、
塑性変形に注意をはらいながら、ここではまずトレンチ
部分の酸化を取り上げ、その酸化形状と応力分布、さら
には所定の形状を得るための設計指針について記述して
みる。例えば、トレンチ角の酸化の「尖り」や「丸め」
現象は、応力が寄与していると言われているが、まだそ
の仕組みが分かっていないのが現状である。
【0412】これらの現象をつぶさに理解するととも
に、上記「尖り」を抑え「丸め」を進行させる的確なプ
ロセス設計を進めることがまず必要である。また、誘起
応力を的確に把握しなければ酸化膜形状も正確には把握
できない。応力の問題は、単に塑性変形上の問題だだけ
でなく、形状把握、さらには、応力場下での不純物分布
把握にも欠かせない項目でもある。また、今後の1G以
降の素子で展開されるであろうSOI基板のLOCOS
工程後の酸化膜厚の予測等は、通常のSi基板と異な
り、酸化進行時にとりわけ応力の影響を受け易いが、そ
の定量的予測はもっと複雑な問題である。
【0413】このような背景において、従来の酸化拡散
モデルに代わり、応力下での酸化剤の挙動把握や、応力
下での酸化現象の理解、応力下での不純物再拡散にまで
作業を広げなければならない。しかも短時間で、高速に
計算されなければならない。本発明者等は、応力の発生
とその再分布の予測、さらには応力下での点欠陥2次元
拡散部分まで計算可能な新しい高速汎用プロセスシミュ
レータ部分を作成し、これを、原子レベルシミュレータ
に付加した。即ち、独自開発による高速2次元酸化/拡
散/応力/変形プロセスシミュレータ部を、原子レベル
材料設計システムの母体に付加した。
【0414】応力のその場観察の実測は、顕微ラマン及
びFT−IRを用いた。これらの計測器のキャリブレー
ションは、あらかじめ行った。SiやSiO2に関する
物性諸量は、工程表が入力された時点で、あらかじめ原
子レベル材料設計システムにより即座に実行され、求め
られる。これが、高速2次元酸化/拡散/応力/変形プ
ロセスシミュレータ部に送られる。
【0415】本発明者らは、まず、この新しい高速2次
元酸化/拡散/応力/変形プロセスシミュレータ部を用
いて突き止めた「尖り」と「丸め」酸化の現象分析を記
述し、更には、適用例として、埋め込み素子分離工程の
基礎検討結果等を示す。Siが酸化されて、SiO2
なると、体積膨張を起こす。このとき、SiO2には粘
弾性があり、その応力の一部は緩和されるものの、依然
SiとSiO2の両方には応力が残留する。
【0416】トレンチコーナーの酸化を例にとり、この
ような状況を調べ、図示したのが図103である。ま
ず、濃度Cの酸化剤がSiO2に到達し、SiO2膜中に
侵入して行くが、この時、酸化剤の拡散D´はSiO2
膜中の応力σの影響を受ける、さらに、Si/SiO2
界面で、酸化剤がSiに反応する時にも基板内の残留蓄
積応力の影響を受ける。
【0417】具体的な効果としては、例えば、SiO2
中に圧縮応力σが存在すると、ネットワークの原子間隔
が狭まり、酸化剤の拡散を抑制することになる。また、
Si基板の面上に圧縮応力σが蓄積すると、ここでも、
Si原子間隔が狭まり、酸化剤の侵入を抑え、結果的に
酸化速度を抑えることになる。
【0418】今、主応力Sのx方向の成分をσxxとし、
また、y方向の成分をσyyとすると、2次元面上の応力
に関しては、下記式(1531)〜(1534)の関係
式が成立する。ここにμは粘性係数、νはポアソン比で
あり、uは変位である。
【0419】
【数141】 ところで、一般に粘性率μの温度依存性は上記式(15
34)で表現できる。ここで、VISC.O、VIS
C.Eは、それぞれ比例定数及び活性化エネルギーを示
す。したがって、上記角応力式にμの温度依存性を加味
すると、温度上昇とともに、各主応力の成分が減少する
ことがわかる。また、Si/SiO2 界面にまで到達し
た酸化剤がSi表面で反応する時の反応定数κs ′は、
同様に上記式(56)のように書き表せる。ここにκs
は、応力が存在しない場合の値である。ここで、σn
は、σn =−(σxxnx 2 +σyyny 2 +2σxynx n
y )であり、στ=−(σxxny 2 +σyynx 2 −2σ
xynx ny )である。VRはパラメータである。この式
から、圧縮応力が増大するとともに、酸化反応定数が減
少することがわかる。上記のSi/SiO2界面で、酸
化剤がSiに反応する時の基板内の残留蓄積応力の影響
がこれである。
【0420】さらに、酸化膜中の酸化剤の拡散係数D′
は、上記式(1531)で表現できる。ここに、pは次
式で表現される。即ちp=−1/2(σxx+σyy)であ
る。またDは応力が存在しない場合の拡散係数である。
この式の意味合いは、以下の様に解釈できる。すなわ
ち、応力σxxやσyyは、引っ張りを正、圧縮を負に取っ
ている。従って、膜内にもし圧縮応力が存在し、物質が
密な状態であれば、拡散係数は小さくなることを示して
いる。上記式は、このことを表現したものである。
【0421】本発明者らが見いだした、この部分のみの
計算手順を図104に示す。まず、Si/SiO2界面
での各節点での酸化剤の濃度を求める。しかる後に、こ
の濃度を基に、1回目の各節点での仮の酸化速度を算出
する。その上で、各節点でのSiO2膜の成長量uを求
める。このuを用いて、上記式(57)から応力Sを求
める。この応力を用いて、式(56)、式(153
1)、さらに式(1534)を用いて酸化速度を再度計
算し直す。そして再び、成長量を求め、応力を算出す
る。そして、十分応力が収束すると、改めて、新しい酸
化時間の進行を始める。即ちニュートン法を用いている
わけである。
【0422】この図104からも分かるように、酸化時
の応力計算は非常に長時間の演算時間を要する。典型的
なトレンチ構造で、初期節点1300点で、900℃の
酸化を試みた。酸化雰囲気はドライ酸素100%とし、
また酸化時間は300分とした。EWS Sun4を用
いて、計算時間を調べてみた。従来技術であれば、概ね
15900分の時間を要する。これはほぼ10日間にも
達し、このままでは、実用には極めて不向きである。
【0423】本発明者等は、ニュートン法部分の計算に
関しては、偏分原理を導入し、前回の収束状況を保存し
ておき、これの1次微分を利用して収束計算を進めた。
この方法により、900℃〜1100℃を例にとり、従
来手法の1000分の1以下に短くすることが出来た。
これにより、本シミュレータは、材料物性値の適正化
と、演算部の高速化を進め、実時間用ツールとして利用
できるところまで到達した。本発明者らはさらに、シミ
ュレータの基本精度として、酸化膜厚の時間依存性値の
精度や、トレンチ角部の形状再現性等を調べてみた。
【0424】図105〜77は、上記シミュレータを用
いて、成長膜厚の温度、時間依存性を調べた結果をグラ
フにまとめたものである。ここには表示しなかったが、
本発明者らは、酸素分圧が10%やさらに1%の時の減
圧挙動についても調べた。その結果、酸化速度が遅いと
きには、酸化膜の膨張に伴う応力が十分緩和しながら進
行するので、基板には応力が蓄積せず、従って、凸部分
にも均一な酸化膜が形成される。本発明者らの慎重な検
討の結果、特に酸化初期に、酸素分圧をさげて、酸化速
度を遅くして、十分緩和させながら徐々に酸素分圧を挙
げて行く方法が有効であることがわかった。
【0425】本発明者らは、本システムを用い、トレン
チ角の「尖り」と「丸め」がなぜ起きるのかを調べた。
既に述べたように、トレンチ角の酸化ではいくつかの要
素を考慮する必要がある。以下に、対応する5個の重要
要素を記す。
【0426】(1)酸化膜粘弾性モデル (2)シリコン弾性モデル (3)酸化剤拡散の応力依存性モデル (4)酸化速度の面方位依存性モデル (5)酸化剤/Si反応のSi表面応力依存性モデル 本発明者らはこれら5個のモデルの内、もっとも支配的
なものは、意外にも第5のモデル即ち、Si/SiO2
界面でのSi基板の応力であることを突き止めた。ここ
で、900℃でのドライ酸化を取り上げることにする。
この時の酸化膜の形状及びSi/SiO2界面の形状を
よくみると、殆どまだ「尖り」形状は認められない。9
00℃ドライ酸化で300分経過後になると、900℃
近傍の酸化温度では、トレンチ角が尖る現象を次のよう
に分析した。
【0427】即ち、100分程度までであれば、さほど
応力は発生しない。しかし、酸化とともに、SiO2
積膨張により、圧縮応力が角部のSiO2膜内とSi先
端部側に集中する。とりわけ、900℃近傍では、粘性
流動効果が小さく、酸化膜の形状緩和効果が小さい。そ
して、Si表面での圧縮応力により、酸化反応が抑止さ
れ、結果として尖るわけである。しかし、1100℃に
なると、粘性流動効果がおおきく、酸化膜の形状緩和効
果がおおきい。そして、Si表面での圧縮応力の蓄積が
なく、結果としてまるくなる。この状況を示したのが図
108である。これらの応力の予測値は、ラマン測定や
IR測定の結果、±10%以内の精度で測定確認でき
た。
【0428】以上の知見を準備しておき、さらに本発明
者らは、上記シミュレーションシステムを、素子分離領
域の面積低減化を目的に、STI(Shallow T
rench Isolarion)の基礎検討の一部に
適用した。ここでは特に、トレンチの上部角の丸め工程
の基礎検討を進めた。酸化温度は1100℃、酸素分圧
は3%で、100分経過後について見た。その結果、以
下のことがわかる。
【0429】(1)はじめは、poly Siの酸化
と、バッファ酸化膜の後退部分の隙間が酸化される。こ
の時、角部は尖り気味である。
【0430】(2)引き続き酸化が進み、隙間が酸化さ
れ、埋め尽くされると、酸化剤の供給が急激に減少し、
側面からの酸化が律速になり、角部は結果的には丸みを
帯びてくる。また主応力は隙間近傍が大きい。
【0431】これに対して、1000℃酸素分圧50%
での酸化時の応力分布について検討した。時間は30分
経過後である。その結果、先の1100℃酸化の応力分
布に比べ、若干応力が大きくなっていることがわかる。
以上の考察から、トレンチ角部の尖りを抑止するには、
Si表面に残留応力を蓄積させないことが第一の条件で
あることがわかる。このための基本的な施策としては、
低温でも十分酸化速度を抑え、酸化膜の変形緩和を十分
起こさせること等が考えられる。
【0432】本発明者らは、さらに塑性変形シミュレー
タとの合体化を進めた。その結果、酸化のみならず、T
EOS膜堆積後の変形等も調べられることがわかった。
【0433】以上の結果をもとに、本発明者らは、特に
すべり面である(111)面上での分解剪断応力に着目
し、2次元シミュレーションシステムでは(110)面
上に投影した値をみた。また熱履歴として、300度/
分であれば良いが、ウエハ内で熱不均一が発生し、部分
的にもし、600度/分にまで達すると、部位によって
は、応力により、転位が発生することも突き止めた。
【0434】以上の予備知見を別途用意しておいて、実
時間でのバイポーラCMOSの製造工程について本発明
を適用した。すなわち、シミュレーションによる逐次予
測と、各工程での必要最小限度のその場測定量との対
比、さらには、各工程の因子のばらつき量と許容量につ
いて検討した。ここではとくに、犠牲酸化の酸化時間ば
らつきと、形状等の許容量とについて実施した例を以下
に示す。
【0435】基本工程は、図109に示す通りである。
さらに引き続き、素子分離用RIEエッチング→酸化→
多結晶Si埋め込み→CMP→多結晶Si部分の犠牲酸
化→の各工程がある。それぞれの工程の因子に変動幅が
ある。この実施例では、多結晶Si部分の犠牲酸化の部
分の特に酸化量に着目した。酸化量の変動幅として、酸
化量零の場合と所定量の2つについて実時間計算結果を
示す。酸化量と、その途上の酸化形状、応力分布さらに
は、この後に続く、LOCOS工程での酸化膜形状及び
その応力の成分特性の変化を示す。
【0436】図110は、1000℃、LOCOS工程
終了時のσxxの分布図、図111は、同じく1000
℃、LOCOS工程終了時のσyyの分布図、また図11
2は、同じく1000℃、LOCOS工程終了時のσxy
の分布図、図113はLOCOS工程終了後、25℃に
温度を下げた時のσxxの分布図、図114は、同じくL
OCOS工程終了後、25℃に温度を下げた時のσyyの
分布図、また図115は、同じくLOCOS工程終了
後、25℃に温度を下げた時のσxyの分布図をそれぞれ
示す。
【0437】これらの図で解るように、CMP後、多結
晶Si膜に対する犠牲酸化が全くない場合は、多結晶S
i膜の表面と窒化膜の表面とが揃っている。この状態で
LOCOS工程を行うと、酸化によって膨張した酸化膜
が窒化膜の横から押すことになる。一方、犠牲酸化を若
干なりとも行った場合は、上手く窒化膜が曲がって行く
ので、応力が緩和して行く。その結果を図116に示し
た。
【0438】また、本発明者らは、本発明者らが提唱し
てきた上記、経験的原子間ポテンシャルを用いた高速分
子動力学シミュレータと汎用プロセスシミュレータとを
合体させたシミュレーションシステムを用いて、酸化工
程や、層間絶縁膜形成工程、さらには、配線工程を経た
後に、MOS素子のSi基板や酸化膜に残留した応力
が、酸化膜の電気的特性に特に大きな影響を及ぼすこと
を突き止めた。
【0439】本発明者らは、上記シミュレータを用い
て、逐次応力生成工程と生成部位をリアルタイムで予測
するとともに、例えば、EEPROM素子部分のトンネ
ル酸化膜が最終的にどのような保持特性等を示すかも出
力させた。さらに、このシステムを用いることにより、
最終的に得られる素子の保持特性の許容幅にたいして、
現在進行中の素子がどのような特性を与えるかも予測で
きる。これが、本発明になるシステムの効力である。さ
らに、工程進捗中に、上手くこのシステムを利用し、最
適な特性を得ることができる様に、引き続く工程を修正
あるいは回避したりすることができる。この実施例を以
下に示す。
【0440】実施例5 この実施例は、本発明をシリコン熱酸化膜中の応力問題
に適用した例を示す。応力問題は、シリコン基板中での
結晶欠陥や転位の発生を抑制するために重要であるが、
それだけではなく、様々な材料において、膜質に影響を
与える要因としても重要である。
【0441】シリコン熱酸化膜には、シリコンが酸化さ
れる際の体積膨張に起因した酸化応力と、昇降温によっ
て生じる熱応力とが生じており、これらが酸化膜質に影
響を与えている。すなわち、原子レベルで言うと、応力
はシリコン酸化膜のネットワーク構造に影響を与え、S
iO2系の単位構造であるテトラヘドロン構造にまで影
響する可能性がある。引いては、電気的な絶縁性にも影
響をもたらす。そこで、本発明者らは、シリコン酸化膜
におけるIRスペクトルや原子レベルでの構造、及び電
気的特性への知見を準備しつつ、酸化膜質を制御しなが
ら酸化工程を進めた。
【0442】最初に、本発明者らは、前段階として、シ
リコン酸化膜のIRスペクトルと原子レベルでの構造及
び電気的特性を調べた。図117に示すように、シリコ
ン(100)基板を用意し、素子分離工程を行って、素
子分離領域により分離された素子領域を作成した。次
に、シリコン基板を酸化装置に導入し、図118に示す
ようにゲート熱酸化膜を形成した。この時、酸化温度を
800℃から1100℃まで種々変化させた。
【0443】所定の膜厚のシリコン酸化膜を形成した
後、基板温度を降温させた。この時、基板温度の降温
は、図119に示すように2段階で行った。第1段階は
基板温度を900℃まで降温する工程であり、基板温度
を急冷させた。第2段階は900℃から室温まで降温さ
せる工程であり、降温速度には特に留意はしなかった。
酸化温度を900℃以下に設定した場合には、第2段階
だけで降温した。
【0444】この降温条件の設定は、別途、酸化膜の粘
性緩和挙動の解析から決定した。一般的に、酸化膜は高
温では柔らかくなり、粘性流動によって応力が緩和しや
すくなると言われている。そこで、本発明者らは、粘性
による応力緩和挙動を詳細に調べた。降温速度による応
力緩和の特性の変化を調べた例を、図120及び図12
1に示した。本例では、酸化温度は1000℃とした。
【0445】図120は典型的な酸化工程を行った後、
それぞれ1000℃で保持した場合、600℃/min
で急冷した場合、6℃/minで徐冷した場合におけ
る、表面に沿った方向の主応力の最大値を調べた結果で
ある。図120に示す結果から、1000℃保持の場合
と徐冷の場合では約1分で応力が緩和し、急冷では粘性
緩和は僅かで、熱応力が重畳されるのみであることが分
かった。
【0446】図121は、同じ降温過程での応力挙動
を、横軸を温度にしてプロットした図である。図121
から、950℃以上の温度領域では応力緩和が顕著であ
り、粘性体もしくは粘弾性体の挙動を示し、950℃以
下では応力緩和が殆ど起こらず、酸化膜は弾性体の挙動
を示すことが分かった。
【0447】このような計算例から推察されるように、
本実施例で行った2段階降温を用いれば、粘性緩和が生
じる温度領域を可能な限り短時間で通過させることによ
り、粘性緩和を抑えることが可能となり、しかも酸化温
度を種々変化させることにより、酸化膜とシリコン基板
の熱膨張係数の差から生じる熱応力を種々設定すること
が出来る。
【0448】このようにして作成した試料を用い、シリ
コン酸化膜直下での基板表面方向の引っ張り応力を顕微
ラマン分光を用いて断面から測定したところ、シリコン
基板中に種々の熱応力が誘起していることが確認され
た。また、本発明者らは、これら基板中の応力誘起測定
値と図122に例示した計算結果とを比較して、酸化膜
中に誘起されている応力値を推察した。尚、本発明者ら
は、これら応力計算を行うに当たって、本発明者らが構
築した分子レベルから汎用レベルまで含めた材料、工程
予測シュミレーションシステムを使用した。その結果、
酸化温度と酸化膜中の圧縮応力との間に、図123に示
す関係が得られた。図123から、酸化温度の上昇とと
もに、酸化膜中の圧縮応力の最大値が直線的に増加する
ことがわかる。
【0449】次の準備として、以上のように形成された
圧縮応力を生じた酸化膜に対して、IRスペクトルを測
定した。また一方、この圧縮応力が付加された場合の酸
化膜構造とIRぺクトルの計算機予測を、古典的分子動
力学法を用いて行った。図124は、圧縮応力付加の有
無による計算IRぺクトルの変化を示している。
【0450】図125は、IRスペクトルの測定値から
Si−O−Si非対称伸縮ピークの応力による波数シフ
トを測定した結果である。実測値と計算値とを比較する
と、圧縮応力でSi−O−Siの非対称伸縮ピークの出
現波数が、応力付加無しの場合よりも低波数側にシフト
する様子が計算IRスペクトルから得られており、実験
と計算は良く一致する事が分かった。そこで、応力を様
々に付加した場合の非晶質構造とIRフルスペクトルを
計算し、準備した。
【0451】さらに次の準備として、応力と電気特性と
の関連を調べた。前述の種々の圧縮応力を生じている酸
化膜をゲート絶縁膜として使用し、その上にゲート電極
を堆積して、Qbdの測定を行った。ゲ一卜電極の堆積
等は、酸化膜の粘性緩和が生じないように、900℃以
下の温度で全ての工程を行つた。そして、前述の分子動
力学計算から予測された酸化膜構造とQbdとの相関を
調べた。
【0452】図126は、種々の応力場での酸化膜のS
i−O原子の動径分布関数を算出した後、第一近接のS
i−O原子間隔が1.8オングストローム以上となって
いる結合の出現頻度とQbdとの相関を求めた結果であ
る。図126から、Si−O原子間隔が1.8A以上と
なっている出現頻度が5%を越えると、Qbdが顕著に
減少しており、Si−0原子間隔と酸化膜絶縁破壊とが
関連していることが分かった。また、IRスペクトルシ
フト量は約15cm-1となつていた。
【0453】酸化膜に圧縮応力が加わった場合の非晶質
ネフトワーク構造の変化は、(1)酸化膜の密度が増加
することから、Si−0結合長は減少するという予想
と、(2)薄膜熱酸化膜のIR測定と中心カモデルか
ら、Si−O−Si結合角度が減少するという予想の2
通りの予想がなされてきた。しかしながら、本発明者ら
は、独自に開発したシミュレータにより、図127に示
したように、酸化膜の密度増加はSi−O結合長の単純
減少にはつながらず、むしろ、結合長分布が広がり、伸
びた結合の割合が増加することを新たに見い出した。加
えて、Qbdとの相関まで見い出した。
【0454】Qbdと伸びたSi−O結合の相関がある
物理的理由については、本発明者らはまだ詳細を解明し
ていないが、その理由のーつを最近報告された論文から
以下のように推察する。引用論文は、「薄膜・表面物理
分科会研究報告、極薄シリコン酸化膜の形成・評価・信
頼性」で金田氏によって報告された「酸化膜中の局所的
変形および不純物による正孔捕獲」p.101−104
である。金田氏が分子軌道法プログラムGAUSSIA
Nを使用して計算したところによると、「伸びたSi−
O−Siがあれば、その中心の酸素への正孔の捕獲が誘
発される。捕獲された正孔は酸素に強く局在するので、
電子に対しては捕獲中心として働く。この結果、酸素の
多くは電子を捕獲して元の状態に戻る。しかし、一部は
電子と正孔の再結合で放出された数eVのエネルギーを
得て、伸びて弱くなっているSi−O結合を切って格子
間へ飛び出し、例えば酸素欠損のような2次的な欠陥を
生成しうる。」と報告されている。
【0455】本発明者らによると、この報告のように、
応力付加により伸びた結合が増加し、酸素欠損欠陥が形
成されることが予想される。さらに、酸素欠損欠陥には
様々な構造が考えられるが、酸素欠損両端のSi原子が
大きく引き離された構造では、電子状態も大きく変化し
て、バンドギャップ内に準位が形成されたり、電気伝導
を生じる電子分布に変化すると推測される。そして酸化
膜の絶縁性が劣化するものと推察される。
【0456】以上のような準備を行うことにより、応力
場中での酸化膜のIRスペクトルが予測でき、さらに、
原子レベルでの構造と電気的特性との相関が得られたこ
とから、本発明者らは、「製造工程の修正」を試みた。
【0457】所定工程まで進んだSiウエハを、酸化膜
形成装置に導入し、1050℃、5%塩酸を含むdry
2雰囲気で、所望の膜厚の酸化膜を形成した。酸化
膜形成後に室温までウエハを急冷し、形成した酸化膜の
IRスペクトルを測定した。測定値より、IRスペクト
ルのSi−O非対称伸縮振動ピーク波数を読みとり、応
力が生じていない場合の同波数からのシフト量を計算し
たところ、17cm-1のシフトが測定された。
【0458】このシフト量から、酸化膜に約0.2GP
aの圧縮応力が酸化膜中に誘起されていることが予測さ
れた。このシフト量をシフト臨界値15cm-1と比較し
た。この熱工程では臨界値を越えていることが検知され
たため、プロセスを停止した。そして、製造工程の修正
として、高温アニールによる応力緩和を行った。
【0459】また、本実施例の効果を比較するため、一
部のウエハには高温アニール工程を付加せず、そのまま
工程を通過させた。応力修正用高温アニールのパラメー
タであるアニール時間とアニール温度、降温速度等は、
予め用意された粘性緩和による緩和時間の解析、及び、
基板中に既に形成されている不純物拡散層の再拡散挙動
を考慮して算出した。この場合は1000℃、30秒の
アニールを行った。
【0460】アニール終了後、再度IRスペクトルを測
定し、シフト量を算出したところ、臨界値を下回ったの
で、そのまま、プロセスを継続した。プロセスを修正し
た場合と、修正しなかった場合とで、Qbdの測定を行っ
た。その結果、修正した場合では、45C/cm2 、修
正しなかった場合で5C/cm2 となり、本発明による
プロセス修正の効果が明確に得られた。
【0461】また、別の構造を有するウエハにおいて、
950℃のdry O2雰囲気で所望の膜厚のゲート酸
化膜を形成した。酸化膜形成後に室温までウエハを急冷
し、形成した酸化膜のIRスペクトルを測定したとこ
ろ、ピークシフト量が16cm-1であった。これは臨界
値を越えていたため、応力緩和のためのアニールの追加
が必要であったが、不純物拡散層の再拡散を考え合わせ
ると、適切なアニール条件の解が存在しなかった。そこ
で、このウエハは工程継続が不可能と判断し、ウエハを
廃棄した。
【0462】
【発明の効果】以上説明したように、本発明によると、
クラスタ化におけるモニタ機能の問題点を解決するだけ
でなく、MOSLSI等の半導体装置を製造する装置、
とりわけ枚葉式あるいはクラスタ化製造装置において、
各装置における限られた量の“その場”測定値や、限ら
れた質的内容の“その場”測定値を最大限に拡張し、テ
ストピース等が無くても、所望の工程通り、あるいは修
正しながら進行するものである。更に、本発明による
と、(i)各プロセス要素の揺らぎや、各要素のずれの
重畳現象を診断し、(ii)所望のシーケンス通り実行さ
れているかのチェックや、(iii)変動因子の発見も可能
である。更にまた、このシステムの特徴は、(iv)プロ
セスの変動を数学的・統計的に、しかも効率よく取り扱
えるようにした点、(v)新しい推論エンジンを付加さ
せた点、及び(vi)新しい逆方向システムを付加した点
にある。
【図面の簡単な説明】
【図1】本発明によるシステムの全体を示す図。
【図2】ポテンシャルの大きさのそれぞれの方向に対す
る傾きが、それぞれの方向の力となることを示す図。
【図3】本発明に用いた新しいabinitio分子動
力学シミュレーションシステムを示す図。
【図4】第一の原子(i)と、第二の原子(j)との作
用に、さらに第三の原子(k)から第二(j)を通して
入ってくる効果を示す図。
【図5】rを粒子間の距離とし、θを粒子角度とし、こ
れらの諸量の主従の関係を示す図。
【図6】従来の計算手法と本発明の高階巡回偏微分法と
の位置づけを示す図。
【図7】分子動力学シミュレータIIと、プロセス変動
算出ユニットIとの連携を説明する図。
【図8】分子動力学シミュレータIIと、プロセス変動
算出ユニットIとの連携を説明する図。
【図9】汎用2・3次元プロセスシミュレータの入力シ
ーケンスの一例を示すを示す図。
【図10】汎用2・3次元プロセスシミュレータの入力
シーケンスの一例を示すを示す図。
【図11】本件発明の実施例による2進木の一部を示す
図。
【図12】イオン化時のクーロン相互作用の模式図を示
す。
【図13】イオン化時のクーロン相互作用の模式図を示
す。
【図14】イオン化時のクーロン相互作用の模式図を示
す。
【図15】簡単な原子構造モデルをを示す。
【図16】アルカリ金属ハライド分子を示す。
【図17】原子価軌道指数ζの表を示す。
【図18】マトリックス(3063)を解いて求めた電荷
を示す。
【図19】マトリックス(3063)を解いて求めた電荷
を示す。
【図20】マトリックス(3063)を解いて求めた電荷
を示す。
【図21】Si-Si間のクーロン相互作用Jの計算結果を示
す。
【図22】Oの電荷が2sにあり、Siが2sあるいは3sにあ
る状態をモデル化した場合のJABの計算結果を示す。
【図23】結晶系及び非晶質系に新たに開発した上記電
荷計算法の手続きをまとめたものを示す。
【図24】本発明による分子動力学シミュレータの実行
結果であるO−Si−O角、Si−O−Si角、及びO
−Si距離を示すデータ図。
【図25】従来のシミュレーションシステムにおける変
分計算を含む各サブルーチンを示す図。
【図26】従来のシミュレーションシステムにおける変
分計算の手順を示す図。
【図27】従来のプロセスと変動プロセスの収束状況を
示す特性図。
【図28】変分の精度を示すデータ図。
【図29】変分の精度を示すデータ図。
【図30】変分の精度を示すデータ図。
【図31】従来のシミュレーションシステムによる不純
物分布とその変動例を示す図。
【図32】従来のシミュレーションシステムによるプロ
セスの変動とその実行組み合わせを示す図。
【図33】従来のシミュレーションシステムによる推論
エンジンの探索の木の一部を示す図。
【図34】本発明者らが開発した新しいabinitioMDシ
ミュレータによるペロブスカイトとSi系結晶の最適化
の図。
【図35】本発明者等による独自開発によるシミュレー
タの構成概念図。
【図36】本発明者等による独自開発によるシミュレー
タの流れ図。
【図37】本発明者等が開発したシミュレータによるP
TOペロブスカイトのC4v対象例の図。
【図38】ペロブスカイトBaTiO3とPbTiO3
構成図。
【図39】本発明者等が新しく作成したシミュレータに
よるペロブスカイトPTOのポテンシャルの出力図。
【図40】本発明者等が開発したシミュレータにおける
SiO2の作成例。
【図41】本発明者等が新しく作成したシミュレータに
よるペロブスカイトPTOのチャージの出力図。
【図42】本発明者等が新しく作成したabinitio/MD
シミュレータによるペロブスカイトPTOの歪みの例。
【図43】本発明者等が新しく作成したabinitio/MD
シミュレータによるペロブスカイトPTOの酸素欠損の
例。
【図44】本発明者等が新しく作成したabinitio/MD
シミュレータによるペロブスカイトPTOのPbの移動
の例。
【図45】本発明者等が新しく作成したabinitio/MD
シミュレータによるペロブスカイトPTOのPbの移動
の例。
【図46】本発明者らが作成したシミュレータのデータ
の一部。
【図47】ソーヤータウワの図。
【図48】本発明の実施例の図。
【図49】本発明に使用されるクラスタ化ツールにおい
て、in−situモニタ測定装置がどのようにして付
加されているかを示す概念図。
【図50】限られた測定量と、他方計算から得られた十
分な情報とを対比させる方法を、モニター上のグラフィ
ックな画像として表示される様子を示す図。
【図51】本発明によるabinitio分子動力学シ
ミュレーションシステムへの入力データの一例としての
温度シーケンスを示す図。
【図52】本発明によるabinitio分子動力学シ
ミュレーションシステムによる原子レベルの実行予測を
示す概念図。
【図53】本発明によるabinitio分子動力学シ
ミュレーションシステムによる原子レベル実行予測を可
視化した概念図。
【図54】本発明によるabinitio分子動力学シ
ミュレーションシステムによるSiO2についての実行
例であるSi−O−Si角の経時変化を示す図。
【図55】図51の入力シーケンスに対する出力例であ
る温度の測定値変化の予測図。
【図56】本発明によるシミュレーションシステムによ
る原子レベル実行例である成膜Si中の配位角の分布
図。
【図57】本発明によるシミュレーションシステムによ
る原子レベル実行例であるSi−Si間距離の分布図。
【図58】本発明によるシミュレーションシステムによ
る原子レベル実行例であるSi−Si間距離の分布図。
【図59】本発明によるシミュレーションシステムによ
る原子レベル実行例であるSi−Si間距離の分布図。
【図60】本発明によるシミュレーションシステムによ
る原子レベル実行例であるSi配位角の分布図。
【図61】本発明によるシミュレーションシステムによ
る原子レベル実行例であるSi−Si間距離の分布図。
【図62】本発明によるシミュレーションシステムによ
る原子レベル実行例。Si−Si間距離の分布図。
【図63】本発明によるシミュレーションシステムによ
る原子レベル実行例であるSi配位角の分布図。
【図64】本発明によるシミュレーションシステムによ
る原子レベル実行例であるSi配位角の分布図。
【図65】本発明によるシミュレーションシステムによ
る原子レベル実行例であるSi配位角の分布図。
【図66】本発明によるシミュレーションシステムによ
る原子レベル実行例であるSi−Si間距離の分布図。
【図67】本発明によるシミュレーションシステムによ
る原子レベル実行例であるSi−Si間距離の分布図。
【図68】本発明によるシミュレーションシステムによ
る原子レベル実行例であるSi−Si間距離の分布図。
【図69】本発明によるシミュレーションシステムによ
る原子レベル実行例であるSi−Si間距離の分布図。
【図70】本発明によって構成された別のシステムの全
体を示す図。
【図71】従来のラマン測定結果を示す特性図。
【図72】従来のラマン測定結果を示す特性図。
【図73】従来のラマン測定結果を示す特性図。
【図74】従来のラマン測定結果を説明するためのモデ
ルを示す図。
【図75】O−Si−O結合角の関数としてのSiO4
四面体の通常モード振動の波長を示す特性図及びストレ
スにより生じた強度の減少を示す特性図。
【図76】SiO2中の欠陥を示すモデル図。
【図77】本発明による分子動力学シミュレータの実行
例を示す特性図。
【図78】本発明による分子動力学シミュレータの実行
例を示す特性図。
【図79】本発明による分子動力学シミュレータの実行
例を示す特性図。
【図80】本発明による分子動力学シミュレータの実行
例である原子の動きを示す図。
【図81】本発明による分子動力学シミュレータの実行
例であるX線スペクトル予測結果を示す特性図。
【図82】本発明による分子動力学シミュレータの実行
例である応力予測結果を示す特性図。
【図83】本発明者らが新しく開発した分子動力学シミ
ュレータによる算出結果を示す特性図。
【図84】本発明による変分手法を用いた分子動力学シ
ミュレータの変動を含む一貫計算を示す図。
【図85】本発明による分子動力学シミュレータによる
探索手法を示す図。
【図86】本発明者らによる非晶質SiO2作成時の温
度設定の様子を示す特性図。
【図87】本発明者らによる計算機中での非晶質SiO
2作成時の圧力変化を示す特性図。
【図88】本発明者らが作成した分子動力学シミュレー
タによる非晶質SiO2作成時の体積変化を示す特性
図。
【図89】本発明者らが新しく開発した分子動力学シミ
ュレータによる非晶質SiO2作成時の体積変化を示す
特性図。
【図90】本発明者らが新しく開発した分子動力学シミ
ュレータによる非晶質SiO2作成時のポテンシャルエ
ネルギー変化を示す特性図。
【図91】本発明者らが新しく開発した分子動力学シミ
ュレータによる非晶質SiO2作成時の計算セル中の2
辺の時間変化を示す特性図。
【図92】本発明者らが新しく開発した分子動力学シミ
ュレータに通常のVerletのアルゴリズムを用い、計算し
たエネルギーの変化を示す特性図。
【図93】本発明者らが新しく開発した分子動力学シミ
ュレータに改良型Verletのアルゴリズムを用い、計算し
たエネルギーの変化を示す特性図。
【図94】本発明者らが新しく開発した分子動力学シミ
ュレータを用い計算したダイポールモーメントの時間変
化を示す特性図。
【図95】本発明者らが新しく開発した分子動力学シミ
ュレータを用い、ダイポールモーメントから更に求めた
パワースペクトルを示す特性図。
【図96】本発明者らが新しく開発した分子動力学シミ
ュレータをα−SiO2 に適用し、計算されたIRスペ
クトルを示す図。
【図97】本発明者らが新しく開発した分子動力学シミ
ュレータをα−SiO2 に適用し、外圧が加わった場合
の予測されたIRスペクトルを示す図。
【図98】本発明者らが新しく開発した分子動力学シミ
ュレータによる非晶質SiO、作成時の体積変化を示す
特性図。
【図99】本発明者らが新しく開発した分子動力学シミ
ュレータに通常のVerletのアルゴリズムを用い、
計算したエネルギーの変化を示す特性図。
【図100】本発明者らによる計算手順を示す図。
【図101】本発明者らによる計算手順を示す図。
【図102】Si(100)における滑り面、分解剪断
応力等を示す図。
【図103】SiO2及びSi/SiO2界面近傍での反
応を示す模式図。
【図104】粘弾性計算の逐次計算を模式的に示す図。
【図105】酸化膜の膜厚の温度、時間依存性を示す特
性図。
【図106】酸化膜の膜厚の温度、時間依存性を示す特
性図。
【図107】酸化膜の膜厚の温度、時間依存性を示す特
性図。
【図108】合わせ込まれた係数を用いた場合のトレン
チ酸化における主応力分布及び最大剪断応力を示す図。
【図109】典型的な熱工程における温度と蓄積応力の
変化を示す特性図。
【図110】トレンチ角におけるσxxの分布を示す計算
予測図。
【図111】トレンチ角におけるσxxの分布を示す計算
予測図。
【図112】トレンチ角におけるσxxの分布を示す計算
予測図。
【図113】トレンチ角におけるσxxの分布を示す計算
予測図。
【図114】トレンチ角におけるσxxの分布を示す計算
予測図。
【図115】トレンチ角におけるσxxの分布を示す計算
予測図。
【図116】最適化された酸化予測の形状を示す図。
【図117】実施例7におけるゲート熱酸化膜の形成工
程を示す断面図。
【図118】実施例7におけるゲート熱酸化膜の形成工
程を示す断面図。
【図119】実施例7におけるゲート熱酸化膜形成後の
昇温工程を示す特性図。
【図120】昇温工程の相違による応力の変化を示す特
性図
【図121】昇温工程の相違による応力の変化を示す特
性図。
【図122】応力分布の計算例を示す図。
【図123】実施例7における酸化温度と膜中の圧縮応
力との関係を示す図。
【図124】応力付加によるIRスペクトルの変化の計
算結果を示す図。
【図125】酸化膜中の圧縮応力とSi−O−Si非対
称伸縮ピークの波数シフトとの関係を示す図。
【図126】第1近傍のSi−O原子間隔が1.8オン
グストローム以上となっている結合の出現頻度とQbdと
の関係を示す図。
【図127】応力付加の有無によるSi−O原子間隔分
布の変化を示す図。
【図128】dRAMの世代とチップコストとの関係を
示す図。
【図129】従来例による半導体装置の製造工程を示す
模式図。
【図130】従来例によるクラスタ化ツールを示す図。
【図131】本発明によるプロセスシミュレーションを
実現する情報処理システムの構成を示す図。
【図132】3次元剛性マトリックスの例を示す図。

Claims (57)

    【特許請求の範囲】
  1. 【請求項1】 電子部品の製造プロセスの製造装置内物
    理現象のシミュレーション方法において、前記製造装置
    内に形成された電子材料の特性に関する製造工程因子を
    入力値とし、この入力値を基に3次元のシミュレーショ
    ンを行ない製造装置内に形成された電子材料の特性の予
    測データを得る方法であって、前記シミュレーション
    は、分子動力学に基づいて行われる工程と、流体モデル
    に基づいて行われる工程とからなり、前記分子動力学に
    基づいて行われるシミュレーション結果として出力され
    る3次元物理量と前記流体モデルに基づいて行われるシ
    ミュレーション結果として出力される3次元物理量の少
    なくとも一方は、他方に転送されそこでのシミュレーシ
    ョンに利用されることを特徴とする電子部品の製造プロ
    セスの設計方法。
  2. 【請求項2】 前記3次元分子動力学シミュレーション
    から前記3次元流体モデルのプロセスシミュレーション
    へ、粘弾性係数、ポアソン比、拡散定数、Cij等の物理
    量が転送されることを特徴とする請求項1に記載の電子
    部品製造プロセスの設計方法。
  3. 【請求項3】 前記3次元流体モデルのプロセスシミュ
    レーションから前記3次元分子動力学シミュレーション
    へ、応力テンソル、温度、不純物濃度、点欠陥濃度等の
    物理量が転送されることを特徴とする請求項1に記載の
    電子部品製造プロセスの設計方法。
  4. 【請求項4】 前記電子部品は、半導体装置であること
    を特徴とする請求項1に記載の電子部品製造プロセスの
    設計方法。
  5. 【請求項5】 更に、前記設定値のゆらぎを入力する段
    階を備え、前記分子動力学シミュレーションは、前記入
    力値とそのゆらぎを基に変分計算を行うことを特徴とす
    る請求項1に記載の電子部品製造プロセスの設計方法。
  6. 【請求項6】 前記分子動力学シミュレーションは、分
    子動力学法と密度汎関数法との組み合わせで行われるこ
    とを特徴とする請求項1に記載の電子部品製造プロセス
    の設計方法。
  7. 【請求項7】 前記分子動力学法は、高階巡回偏微分法
    に基づいていることを特徴とする請求項6に記載の電子
    部品製造プロセスの設計方法。
  8. 【請求項8】 前記高階巡回偏微分法は、3体問題を扱
    っていることを特徴とする請求項7に記載の電子部品製
    造プロセスの設計方法。
  9. 【請求項9】 前記分子動力学シミュレーションは、経
    験的ポテンシャルを用いて行われることを特徴とする請
    求項1に記載の電子部品製造プロセスの設計方法。
  10. 【請求項10】 前記経験的ポテンシャルはATTAポ
    テンシャルであることを特徴とする請求項9に記載の電
    子部品製造プロセスの設計方法。
  11. 【請求項11】 前記経験的ポテンシャルはBKSポテ
    ンシャルであることを特徴とする請求項9に記載の電子
    部品製造プロセスの設計方法。
  12. 【請求項12】 電子部品の製造工程において、製造装
    置内に形成された電子材料の特性を測定して実観測デー
    タを得る工程と、前記製造装置内に形成された電子材料
    の特性に関する製造工程因子を入力値とし、この入力値
    を基に3次元のシミュレーションを行ない前記製造装置
    内に形成された電子材料の特性の予測データを得る工程
    と、前記予測データと実観測データとを比較検定する工
    程と、前記比較検定により、前記製造工程因子の設定値
    と、前記実観測データから推測される前記製造工程因子
    との間に有意差が認められた場合、前記製造工程因子を
    修正処理する工程とを備え、前記シミュレーションは、
    分子動力学に基づいて行われる工程と、流体モデルに基
    づいて行われる工程とからなり、前記分子動力学に基づ
    いて行われるシミュレーション結果として出力される3
    次元物理量と前記流体モデルに基づいて行われるシミュ
    レーション結果として出力される3次元物理量の少なく
    とも一方は、他方に転送されそこでのシミュレーション
    に利用されることを特徴とする電子部品の製造方法。
  13. 【請求項13】 前記実観測データから推測される前記
    製造工程因子との間の有意差が許容範囲を越えていて、
    修正処理が不能である場合には、不良品として処理する
    ことを特徴とする請求項12に記載の電子部品の製造方
    法。
  14. 【請求項14】 前記製造装置はクラスタ化しているこ
    とを特徴とする請求項12に記載の電子部品の製造方
    法。
  15. 【請求項15】 前記クラスタ化した製造装置は枚葉式
    であることを特徴とする請求項14に記載の電子部品の
    製造方法。
  16. 【請求項16】 前記3次元分子動力学シミュレーショ
    ンから前記3次元流体モデルのプロセスシミュレーショ
    ンへ、粘弾性係数、ポアソン比、拡散定数、Cij等の物
    理量が転送されることを特徴とする請求項12に記載の
    電子部品の製造方法。
  17. 【請求項17】 前記3次元流体モデルのプロセスシミ
    ュレーションから前記3次元分子動力学シミュレーショ
    ンへ、応力テンソル、温度、不純物濃度、点欠陥濃度等
    の物理量が転送されることを特徴とする請求項12に記
    載の電子部品の製造方法。
  18. 【請求項18】 前記特性は、半導体基板上に形成され
    た膜、半導体基板表面、及び半導体基板内部の少なくと
    も1つに関する特性であることを特徴とする請求項12
    に記載の電子部品の製造方法。
  19. 【請求項19】 前記予測データと実観測データとの比
    較検定は逐次実時間で行われることを特徴とする請求項
    12に記載の電子部品の製造方法。
  20. 【請求項20】 前記製造工程因子を修正処理する工程
    は逐次実時間で行われることを特徴とする請求項12に
    記載の電子部品の製造方法。
  21. 【請求項21】 前記分子動力学シミュレーションは、
    分子動力学法と密度汎関数法との組み合わせで行われる
    ことを特徴とする請求項12に記載の電子部品の製造方
    法。
  22. 【請求項22】 前記分子動力学法は、高階巡回偏微分
    法に基づいていることを特徴とする請求項21に記載の
    電子部品の製造方法。
  23. 【請求項23】 前記高階巡回偏微分法は、3体問題を
    扱っていることを特徴とする請求項22に記載の電子部
    品の製造方法。
  24. 【請求項24】 前記分子動力学シミュレーションは、
    経験的ポテンシャルを用いて行われることを特徴とする
    請求項12に記載の電子部品の製造方法。
  25. 【請求項25】 前記経験的ポテンシャルはATTAポ
    テンシャルであるすることを特徴とする請求項24に記
    載の電子部品の製造方法。
  26. 【請求項26】 前記経験的ポテンシャルはBKSポテ
    ンシャルであるすることを特徴とする請求項24に記載
    の電子部品の製造方法。
  27. 【請求項27】 前記分子動力学シミュレーションで得
    られた前記複数の特性の少なくとも1つの予測データを
    用いて、更に流体モデルのプロセスシミュレーションが
    行われることを特徴とする請求項12に記載の電子部品
    の製造方法。
  28. 【請求項28】 前記分子動力学シミュレーションと流
    体モデルのプロセスシミュレーションは、共に3次元シ
    ミュレーションであることを特徴とする請求項27に記
    載の電子部品の製造方法。
  29. 【請求項29】 前記3次元分子動力学シミュレーショ
    ンから前記3次元流体モデルのプロセスシミュレーショ
    ンへ、粘弾性係数、ポアソン比、拡散定数、Cij等の物
    理量が転送されることを特徴とする請求項27に記載の
    電子部品の製造方法。
  30. 【請求項30】 電子部品の製造装置と、前記製造装置
    内に形成された前記電子部品の電子材料の特性を測定す
    る測定装置と、前記電子部品の製造工程因子の設定値を
    入力値とし、この入力値を基にシミュレーションを行な
    い前記特性の予測データを得るコンピュータシステム
    と、前記測定装置及び前記製造装置と前記コンピュータ
    システムを接続する手段とを備え、前記コンピュータシ
    ステムは、前記予測データと実観測データとを比較検定
    を行ない、前記比較検定により、前記製造工程因子の設
    定値と、前記実観測データから推測される前記製造工程
    因子との間に有意差が認められた場合、前記製造装置の
    制御手段を介して製造工程因子を修正処理することを特
    徴とする電子部品の製造システム。
  31. 【請求項31】 前記製造装置はクラスタ化しているこ
    とを特徴とする請求項30に記載の電子部品の製造シス
    テム。
  32. 【請求項32】 前記クラスタ化した製造装置は枚葉式
    であることを特徴とする請求項31に記載の電子部品の
    製造システム。
  33. 【請求項33】 前記シミュレーションは、分子動力学
    法と密度汎関数法との組み合わせによる分子動力学シミ
    ュレーションに基づいて行われることを特徴とする請求
    項30に記載の電子部品の製造システム。
  34. 【請求項34】 前記特性は、半導体基板上に形成され
    た膜、半導体基板表面、及び半導体基板内部の少なくと
    も1つに関する特性であることを特徴とする請求項30
    に記載の電子部品の製造システム。
  35. 【請求項35】 前記予測データと実観測データとの比
    較検定は逐次実時間で行われることを特徴とする請求項
    30に記載の電子部品の製造システム。
  36. 【請求項36】 前記製造工程因子を修正処理する工程
    は逐次実時間で行われることを特徴とする請求項30に
    記載の電子部品の製造システム。
  37. 【請求項37】 更に、前記入力値として前記設定値の
    ゆらぎが入力され、前記分子動力学シミュレーション
    は、前記入力値とそのゆらぎを基に変分計算を行うこと
    を特徴とする請求項30に記載の電子部品の製造システ
    ム。
  38. 【請求項38】 前記分子動力学法は、高階巡回偏微分
    法に基づいていることを特徴とする請求項30に記載の
    電子部品の製造システム。
  39. 【請求項39】 前記高階巡回偏微分法は、3体問題を
    扱っていることを特徴とする請求項38に記載の電子部
    品の製造システム。
  40. 【請求項40】 前記分子動力学シミュレーションは、
    経験的ポテンシャルを用いた方法を併用することを特徴
    とする請求項30に記載の電子部品の製造システム。
  41. 【請求項41】 前記経験的ポテンシャルはATTAポ
    テンシャルであるすることを特徴とする請求項40に記
    載の電子部品の製造システム。
  42. 【請求項42】 前記経験的ポテンシャルはBKSポテ
    ンシャルであるすることを特徴とする請求項40に記載
    の電子部品の製造システム。
  43. 【請求項43】 前記分子動力学シミュレーションで得
    られた前記特性の少なくとも1つの予測データを用い
    て、更に流体モデルのプロセスシミュレーションが行わ
    れることを特徴とする請求項30に記載の電子部品の製
    造システム。
  44. 【請求項44】 前記分子動力学シミュレーションと流
    体モデルのプロセスシミュレーションは、共に3次元シ
    ミュレーションであることを特徴とする請求項43に記
    載の電子部品の製造システム。
  45. 【請求項45】 前記3次元分子動力学シミュレーショ
    ンから前記3次元流体モデルのプロセスシミュレーショ
    ンへ、粘弾性係数、ポアソン比、拡散定数、Cij等の物
    理量が転送されることを特徴とする請求項43に記載の
    電子部品の製造システム。
  46. 【請求項46】 電子部品の製造プロセスの製造装置内
    物理現象のシミュレーションを行うプログラムであっ
    て、前記製造装置内に形成された電子材料の特性に関す
    る製造工程因子の設定値を入力する手順と、この入力値
    を基に分子動力学法と密度汎関数法との組み合わせによ
    る分子動力学シミュレーションを行い、前記特性の予測
    データを得る手順とを実行させるプログラムを記録した
    コンピュータ読み取り可能な記録媒体。
  47. 【請求項47】 前記分子動力学シミュレーションは、
    分子動力学法と密度汎関数法との組み合わせによること
    を特徴とする請求項46に記載の記録媒体。
  48. 【請求項48】 前記分子動力学法は、高階巡回偏微分
    法に基づいていることを特徴とする請求項47に記載の
    記録媒体。
  49. 【請求項49】 前記高階巡回偏微分法は、3体問題を
    扱っていることを特徴とする請求項48に記載の記録媒
    体。
  50. 【請求項50】 前記電子部品は、半導体装置であるこ
    とを特徴とする請求項47に記載の記録媒体。
  51. 【請求項51】 更に、前記設定値のゆらぎを入力する
    手順を備え、前記分子動力学シミュレーションは、前記
    入力値とそのゆらぎを基に変分計算を行うことを特徴と
    する請求項47に記載の記録媒体。
  52. 【請求項52】 前記分子動力学シミュレーションは、
    経験的ポテンシャルを用いた方法を併用することを特徴
    とする請求項47に記載の記録媒体。
  53. 【請求項53】 前記経験的ポテンシャルはATTAポ
    テンシャルであるすることを特徴とする請求項52に記
    載の記録媒体。
  54. 【請求項54】 前記経験的ポテンシャルはBKSポテ
    ンシャルであるすることを特徴とする請求項52に記載
    の記録媒体。
  55. 【請求項55】 前記分子動力学シミュレーションで得
    られた前記特性の予測データを用いて、更に流体モデル
    のプロセスシミュレーションが行われることを特徴とす
    る請求項47に記載の記録媒体。
  56. 【請求項56】 前記分子動力学シミュレーションと流
    体モデルのプロセスシミュレーションは、共に3次元シ
    ミュレーションであることを特徴とする請求項55に記
    載の記録媒体。
  57. 【請求項57】 前記3次元分子動力学シミュレーショ
    ンから前記3次元流体モデルのプロセスシミュレーショ
    ンへ、粘弾性係数、ポアソン比、拡散定数、Cij等の物
    理量が転送されることを特徴とする請求項55に記載の
    記録媒体。
JP34644097A 1997-12-16 1997-12-16 電子部品の製造方法、製造システム、設計方法、及び記録媒体 Withdrawn JPH11176906A (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP34644097A JPH11176906A (ja) 1997-12-16 1997-12-16 電子部品の製造方法、製造システム、設計方法、及び記録媒体

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP34644097A JPH11176906A (ja) 1997-12-16 1997-12-16 電子部品の製造方法、製造システム、設計方法、及び記録媒体

Publications (1)

Publication Number Publication Date
JPH11176906A true JPH11176906A (ja) 1999-07-02

Family

ID=18383447

Family Applications (1)

Application Number Title Priority Date Filing Date
JP34644097A Withdrawn JPH11176906A (ja) 1997-12-16 1997-12-16 電子部品の製造方法、製造システム、設計方法、及び記録媒体

Country Status (1)

Country Link
JP (1) JPH11176906A (ja)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007507889A (ja) * 2003-09-30 2007-03-29 東京エレクトロン株式会社 半導体処理ツールによって実行されるプロセスを分析する第1の原理シミュレーションを使用するシステム及び方法。
JP2007507886A (ja) * 2003-09-30 2007-03-29 東京エレクトロン株式会社 半導体製造プロセスを制御する第1の原理シミュレーションを用いたシステム及び方法。
JP2007507890A (ja) * 2003-09-30 2007-03-29 東京エレクトロン株式会社 半導体製造プロセスにおいて第1の原理シミュレーションを用いたシステム及び方法。
JP2007507888A (ja) * 2003-09-30 2007-03-29 東京エレクトロン株式会社 半導体製造プロセスを制御するために第1の原理シミュレーションを用いたシステム及び方法。
WO2011036952A1 (ja) * 2009-09-28 2011-03-31 株式会社日立製作所 電子状態計算システム及びプログラム
US8032348B2 (en) 2003-09-30 2011-10-04 Tokyo Electron Limited System and method for using first-principles simulation to facilitate a semiconductor manufacturing process
KR20190115895A (ko) * 2018-04-04 2019-10-14 세메스 주식회사 확산계수 시뮬레이션 장치 및 확산계수 시뮬레이션 방법
JP2020149617A (ja) * 2019-03-15 2020-09-17 学校法人慶應義塾 分子動力学データ解析装置及びプログラム
CN114840800A (zh) * 2022-04-11 2022-08-02 武汉大学 一种确定晶圆键合过程中最佳键合力的理论方法

Cited By (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8296687B2 (en) 2003-09-30 2012-10-23 Tokyo Electron Limited System and method for using first-principles simulation to analyze a process performed by a semiconductor processing tool
JP2007507890A (ja) * 2003-09-30 2007-03-29 東京エレクトロン株式会社 半導体製造プロセスにおいて第1の原理シミュレーションを用いたシステム及び方法。
US8050900B2 (en) 2003-09-30 2011-11-01 Tokyo Electron Limited System and method for using first-principles simulation to provide virtual sensors that facilitate a semiconductor manufacturing process
US8073667B2 (en) 2003-09-30 2011-12-06 Tokyo Electron Limited System and method for using first-principles simulation to control a semiconductor manufacturing process
JP4795957B2 (ja) * 2003-09-30 2011-10-19 東京エレクトロン株式会社 半導体製造プロセスを制御する第1の原理シミュレーションを用いたシステム及び方法。
US8014991B2 (en) 2003-09-30 2011-09-06 Tokyo Electron Limited System and method for using first-principles simulation to characterize a semiconductor manufacturing process
US8032348B2 (en) 2003-09-30 2011-10-04 Tokyo Electron Limited System and method for using first-principles simulation to facilitate a semiconductor manufacturing process
US8036869B2 (en) 2003-09-30 2011-10-11 Tokyo Electron Limited System and method for using first-principles simulation to control a semiconductor manufacturing process via a simulation result or a derived empirical model
JP2007507889A (ja) * 2003-09-30 2007-03-29 東京エレクトロン株式会社 半導体処理ツールによって実行されるプロセスを分析する第1の原理シミュレーションを使用するシステム及び方法。
JP2007507886A (ja) * 2003-09-30 2007-03-29 東京エレクトロン株式会社 半導体製造プロセスを制御する第1の原理シミュレーションを用いたシステム及び方法。
JP2007507888A (ja) * 2003-09-30 2007-03-29 東京エレクトロン株式会社 半導体製造プロセスを制御するために第1の原理シミュレーションを用いたシステム及び方法。
KR101094620B1 (ko) 2003-09-30 2011-12-15 도쿄엘렉트론가부시키가이샤 반도체 프로세싱 도구에 의해 수행되는 프로세스를 용이하게 하는 방법 및 시스템, 시스템, 및 컴퓨터 판독가능한 매체
JP5430665B2 (ja) * 2009-09-28 2014-03-05 株式会社日立製作所 電子状態計算システム及びプログラム
WO2011036952A1 (ja) * 2009-09-28 2011-03-31 株式会社日立製作所 電子状態計算システム及びプログラム
KR20190115895A (ko) * 2018-04-04 2019-10-14 세메스 주식회사 확산계수 시뮬레이션 장치 및 확산계수 시뮬레이션 방법
JP2020149617A (ja) * 2019-03-15 2020-09-17 学校法人慶應義塾 分子動力学データ解析装置及びプログラム
CN114840800B (zh) * 2022-04-11 2024-07-16 武汉大学 一种确定晶圆键合过程中最佳键合力的理论方法
CN114840800A (zh) * 2022-04-11 2022-08-02 武汉大学 一种确定晶圆键合过程中最佳键合力的理论方法

Similar Documents

Publication Publication Date Title
US6185472B1 (en) Semiconductor device manufacturing method, manufacturing apparatus, simulation method and simulator
Brehm et al. Tunable quadruple-well ferroelectric van der Waals crystals
Morozovska et al. Nanoscale electromechanics of paraelectric materials with mobile charges: Size effects and nonlinearity of electromechanical response of SrTiO 3 films
Rabe et al. First-principles studies of ferroelectric oxides
Paulauskas et al. Atomic scale study of polar Lomer–Cottrell and Hirth lock dislocation cores in CdTe
Braithwaite et al. A theoretical study of the energetics and IR frequencies of hydroxyl defects in forsterite
JPH11176906A (ja) 電子部品の製造方法、製造システム、設計方法、及び記録媒体
JP4516030B2 (ja) 半導体装置の製造プロセスの設計方法、半導体装置の製造方法、及び半導体装置の製造システム
Deluca et al. Advantages and developments of Raman spectroscopy for electroceramics
Tharpe et al. Resolving mechanical properties and morphology evolution of free‐standing ferroelectric Hf0. 5Zr0. 5O2
Chae et al. Local epitaxial templating effects in ferroelectric and antiferroelectric ZrO2
JP3866348B2 (ja) 半導体装置の製造方法、製造装置、シミュレーション方法、及びシミュレータ
Jaszewski et al. Infrared Signatures for Phase Identification in Hafnium Oxide Thin Films
JP2006013525A (ja) 半導体装置の製造方法、製造装置、シミュレーション方法、及びシミュレータ
Lin et al. Strain effects on domain structures in ferroelectric thin films from phase‐field simulations
Ahart et al. Brillouin scattering and molecular dynamics study of the elastic properties of Pb (Mg 1∕ 3 Nb 2∕ 3) O 3
Vamvakas et al. Comparison of FTIR transmission spectra of thermally and LPCVD SiO2 films grown by TEOS pyrolysis
Ouyang et al. Quantum-Accurate Modeling of Ferroelectric Phase Transition in Perovskites from Message-Passing Neural Networks
Alfredsson et al. Crystal morphology and surface structures of orthorhombic MgSiO 3 perovskite
Cervasio et al. IR Spectroscopy of PbZr1–x Ti x O3 Material: A Complementary Experimental/Ab Initio Calculations Approach of a Solid Solution
Zatryb et al. Raman scattering from confined acoustic phonons of silicon nanocrystals in silicon oxide matrix
De La Pierre et al. Vibrational spectroscopy of minerals through ab initio methods
Belboukhari et al. Efficient exploration of electronic and dielectric properties using advanced
Varjas et al. Computation of topological invariants of disordered materials using the kernel polynomial method
Vaideeswaran In Search of Ferroelectricity in Antiferroelectric Lead Zirconate

Legal Events

Date Code Title Description
A761 Written withdrawal of application

Free format text: JAPANESE INTERMEDIATE CODE: A761

Effective date: 20051226