JP4192704B2 - A simulation method to maximize the defect-free region of a single crystal - Google Patents
A simulation method to maximize the defect-free region of a single crystal Download PDFInfo
- Publication number
- JP4192704B2 JP4192704B2 JP2003194779A JP2003194779A JP4192704B2 JP 4192704 B2 JP4192704 B2 JP 4192704B2 JP 2003194779 A JP2003194779 A JP 2003194779A JP 2003194779 A JP2003194779 A JP 2003194779A JP 4192704 B2 JP4192704 B2 JP 4192704B2
- Authority
- JP
- Japan
- Prior art keywords
- single crystal
- temperature
- distribution
- mesh
- defect
- 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.)
- Expired - Lifetime
Links
Images
Landscapes
- Crystals, And After-Treatments Of Crystals (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、チョクラルスキー（以下、ＣＺという。）法にて引上げられるシリコン等の単結晶の無欠陥領域をコンピュータシミュレーションにより最大化する方法に関するものである。
【０００２】
【従来の技術】
近年の半導体集積回路の超微細化にともないデバイスの歩留まりを低下させる要因として、結晶に起因したパーティクル（Crystal Originated Particle、以下、ＣＯＰという。）や、酸化誘起積層欠陥（Oxidation Induced Stacking Fault、以下、ＯＩＳＦという。）の核となる酸素析出物の微小欠陥や、或いは侵入型転位（Interstitial-type Large Dislocation、以下、Ｌ／Ｄという。）の存在が挙げられている。
【０００３】
ＣＯＰは、鏡面研磨されたシリコンウェーハをアンモニアと過酸化水素の混合液で洗浄すると、ウェーハ表面に出現するピットである。このウェーハをパーティクルカウンタで測定すると、このピットがパーティクルとして検出される。このピットは結晶に起因する。ＣＯＰは電気的特性、例えば酸化膜の経時絶縁破壊特性（Time Dependent dielectric Breakdown、ＴＤＤＢ）、酸化膜耐圧特性（Time Zero Dielectric Breakdown、ＴＺＤＢ）等を劣化させる原因となる。またＣＯＰがウェーハ表面に存在するとデバイスの配線工程において段差を生じ、断線の原因となり得る。そして素子分離部分においてもリーク等の原因となり、製品の歩留りを低くする。
【０００４】
ＯＩＳＦは、結晶成長時に微小な酸素析出核がシリコン単結晶中に導入され、半導体デバイスを製造する際の熱酸化工程等で顕在化する欠陥である。このＯＩＳＦは、デバイスのリーク電流を増加させる等の不良原因になる。Ｌ／Ｄは、転位クラスタとも呼ばれたり、或いはこの欠陥を生じたシリコンウェーハをフッ酸を主成分とする選択エッチング液に浸漬するとピットを生じることから転位ピットとも呼ばれる。このＬ／Ｄも、電気的特性、例えばリーク特性、アイソレーション特性等を劣化させる原因となる。以上のことから、半導体集積回路を製造するために用いられるシリコンウェーハからＣＯＰ、ＯＩＳＦ及びＬ／Ｄを減少させることが必要となっている。
【０００５】
このＣＯＰ、ＯＩＳＦ及びＬ／Ｄを有しない無欠陥の単結晶インゴット及びこの単結晶インゴットからスライスされたシリコンウェーハが開示されている（例えば、特許文献１参照）。この無欠陥の単結晶インゴットは、単結晶インゴット内での空孔型点欠陥の凝集体及び格子間シリコン型点欠陥の凝集体がそれぞれ存在しないパーフェクト領域を［Ｐ］とするとき、パーフェクト領域［Ｐ］からなる単結晶インゴットである。パーフェクト領域［Ｐ］は、単結晶インゴット内で空孔型点欠陥が優勢であって凝集した点欠陥を有する領域［Ｖ］と、格子間シリコン型点欠陥が優勢であって凝集した点欠陥を有する領域［Ｉ］との間に介在する。
【０００６】
【特許文献１】
米国特許第６，０４５，６１０号明細書に対応する特開平１１−１３９３号公報
【０００７】
【発明が解決しようとする課題】
しかし、上記従来の特開平１１−１３９３号公報に示された無欠陥の単結晶インゴットでは、単結晶インゴットの全長にわたって空孔型点欠陥の凝集体及び格子間シリコン型点欠陥の凝集体が存在しない引上げ条件を設定するために、引上げ機の種類や単結晶インゴットの太さ毎に実際に単結晶インゴットを引上げてそのライフタイムを測定するという実験を数多く繰返さなければならず、最適な引上げ条件を設定するのに極めて多くの時間と労力を要する不具合があった。
また、上記従来の無欠陥の単結晶インゴットでは、点欠陥濃度を求めるために用いているＶ／Ｇ理論、即ち単結晶の引上げ速度をＶ（ｍｍ／分）とし、単結晶及び融液の固液界面近傍の温度勾配をＧ（℃／ｍｍ）とするときに、Ｖ／Ｇ（ｍｍ2／分・℃）を制御するボロンコフの理論を簡略化した一次元モデル（単結晶の引上げ方向だけのモデル）であるため、単結晶の半径方向への点欠陥の拡散（外方拡散）が考慮されておらず、単結晶の半径方向への点欠陥の濃度勾配を求めることができない問題点があった。
更に、上記従来の無欠陥の単結晶インゴットでは、Ｖ／Ｇ理論にて点欠陥の濃度を求める場合に、上記Ｇを高精度で求めるために融液の対流をモデル化して計算する必要があるにも拘らず、このような考慮が全くなされていないため、単結晶内の点欠陥の濃度分布を精度良く求めることができない問題点もあった。
【０００８】
本発明の目的は、実際に単結晶を引上げるという実験の回数を少なくすることができ、これにより比較的少ない時間と労力で、単結晶の引上げ方向及び半径方向に無欠陥領域が最も拡大する単結晶製造条件を計算により正確に求めることができる、シミュレーション方法を提供することにある。
【０００９】
【課題を解決するための手段】
請求項１に係る発明は、図１〜図７に示すように、引上げ機１１により単結晶１４を引上げるときの単結晶製造条件としてフィードバックされる変数を一定の間隔で変量させたパラメータ群Ｐ1，Ｐ2，…，ＰNを任意に定義する第１ステップと、パラメータ群Ｐ1，Ｐ2，…，ＰNのうちパラメータＰ1を単結晶製造条件として単結晶１４の融液からの引上げ開始時から単結晶１４の冷却完了時までの引上げ機１１のホットゾーンをメッシュ構造でモデル化する第２ステップと、ホットゾーンの各部材毎にメッシュをまとめかつこのまとめられたメッシュに対する各部材の物性値とともにパラメータ群Ｐ1，Ｐ2，…，ＰNのうちのパラメータＰ1をコンピュータに入力する第３ステップと、各部材の表面温度分布をヒータの発熱量及び各部材の輻射率に基づいて求める第４ステップと、各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式を解くことにより各部材の内部温度分布を求めた後に融液１２が乱流であると仮定して得られた乱流モデル式及びナビエ・ストークスの方程式を連結して解くことにより対流を考慮した融液１２の内部温度分布を更に求める第５ステップと、単結晶１４及び融液１２の固液界面形状を単結晶の三重点Ｓを含む等温線に合せて求める第６ステップと、第４ステップから第６ステップを三重点Ｓが単結晶１４の融点になるまで繰返し引上げ機１１内の温度分布を計算して単結晶１４のメッシュの座標及び温度を求めこれらのデータをそれぞれコンピュータに入力する第７ステップと、単結晶１４の引上げ長及び引上げ高さを段階的に変えて第２ステップから第７ステップまでを繰返し引上げ機１１内の温度分布を計算して単結晶１４のメッシュの座標及び温度を求めこれらのデータをそれぞれコンピュータに入力する第８ステップと、単結晶１４の融液１２からの引上げ開始時から単結晶１４の冷却完了時までの時間を所定の間隔毎に区切りこの区切られた時間間隔毎に第８ステップで求めた単結晶１４のメッシュの座標及び温度のデータから単結晶１４の引上げ長及び引上げ高さと単結晶１４内の温度分布とを求める第９ステップと、単結晶１４内の空孔及び格子間原子の拡散係数及び境界条件に基づいて拡散方程式を解くことにより所定の時間間隔の経過した後の空孔及び格子間原子の濃度分布を求める第１０ステップと、空孔の濃度分布に基づいてボイドの形成開始温度を求める第１１ステップと、空孔の濃度分布に基づいて高温酸素析出物の形成開始温度を求める第１２ステップと、単結晶１４内のそれぞれのメッシュの格子点における温度が次第に低下してボイドの形成開始温度になったときのボイドの濃度を求める第１３ステップと、単結晶１４内のそれぞれのメッシュの格子点における温度がボイドの形成開始温度より低いときのボイドの半径を求める第１４ステップと、単結晶１４内のそれぞれのメッシュの格子点における温度が次第に低下して高温酸素析出物の形成開始温度になったときの高温酸素析出物の濃度を求める第１５ステップと、単結晶１４内のそれぞれのメッシュの格子点における温度が高温酸素析出物の形成開始温度より低いときの高温酸素析出物の半径を求める第１６ステップと、第９ステップから第１６ステップを単結晶１４の冷却が完了するまで繰返す第１７ステップと、高温酸素析出物の存在するＰバンドとボイド及び高温酸素析出物のいずれも存在しない無欠陥領域とを区画しかつ無欠陥領域がＰバンドより単結晶ボトム側に位置するように区画する第１等濃度線ＨＣ１Xを求める第１８ステップと、高濃度酸素析出物が存在するＢバンドと無欠陥領域とを区画しかつ無欠陥領域がＢバンドより単結晶トップ側に位置するように区画する第１分布線ＢＣ１Xを求める第１９ステップと、第１等濃度線ＨＣ１Xの変曲点の最大値と第１分布線ＢＣ１Xの変曲点の最小値との差ΔＺ1を求める第２０ステップと、単結晶製造条件のパラメータをＰ2からＰNに順次変え第２ステップから第２０ステップを繰返して第１等濃度線ＨＣ１Xの変曲点の最大値と第１分布線ＢＣ１Xの変曲点の最小値との差ΔＺ1を求めた後にこの差ΔＺ1が最も大きくなる単結晶製造条件を求める第２１ステップとを含むコンピュータを用いて単結晶の無欠陥領域を最大化するシミュレーション方法である。
【００１０】
この請求項１に記載された単結晶の無欠陥領域を最大化するシミュレーション方法では、パラメータＰ1の単結晶製造条件で融液１２の対流を考慮して融液１２から成長する単結晶１４内の温度分布を求めるだけでなく、更に冷却過程における単結晶１４内の温度分布までも求めることによって、即ち融液１２から切離された単結晶１４の冷却過程における単結晶１４の徐冷及び急冷の効果を考慮して解析することによって、単結晶１４内のボイドの濃度分布及びサイズ分布を予測するとともに、単結晶１４内の高温酸素析出物の濃度分布及びサイズ分布も予測する。次いで単結晶１４内の第１等濃度線ＨＣ１Xと第１分布線ＢＣ１Xを計算により求めた後に、第１等濃度線ＨＣ１Xの変曲点の最大値と第１分布線ＢＣ１Xの変曲点の最小値との差ΔＺ1を計算により求める。
次に上記単結晶製造条件のパラメータをＰ2から順にＰNまで変え、上記と同様にして第１等濃度線ＨＣ１Xの変曲点の最大値と第１分布線ＢＣ１Xの変曲点の最小値との差ΔＺ1を計算により求める。更に上記第１等濃度線ＨＣ１Xの変曲点の最大値と第１分布線ＢＣ１Xの変曲点の最小値との差ΔＺ1が最も大きくなる単結晶製造条件を求める。これにより単結晶１４の引上げ方向及び半径方向に無欠陥領域が最も拡大する単結晶製造条件を計算により正確に求めることができる。
【００１１】
なお、本明細書において「高温酸素析出物」とは、ボイド欠陥の内側に析出する内壁酸化膜を球形に近似して扱ったものをいう。また酸素析出物を高温と限定したのは、ボイド欠陥より低温で析出するＢＭＤ（Bulk Micro Defect）等の析出核となるものと区別するためである。また、本明細書において「無欠陥領域」には、空孔が優勢な真性点欠陥を含む無欠陥領域及び格子間原子が優勢な真性欠陥を含む無欠陥領域が存在する。
【００１２】
請求項２に係る発明は、図６及び図８〜図１２に示すように、請求項１に記載の第１ステップから第１７ステップと、高温酸素析出物の存在するＰバンドとボイド及び高温酸素析出物のいずれも存在しない無欠陥領域とを区画しかつ無欠陥領域がＰバンドより単結晶トップ側に位置するように区画する第２等濃度線ＨＣ２Xを求める第１８ステップと、高濃度酸素析出物が存在するＢバンドと無欠陥領域とを区画しかつ無欠陥領域がＢバンドより単結晶ボトム側に位置するように区画する第２分布線ＢＣ２Xを求める第１９ステップと、第２等濃度線ＨＣ２Xの変曲点の最小値と第２分布線ＢＣ２Xの変曲点の最大値との差ΔＺ2を求める第２０ステップと、単結晶製造条件のパラメータをＰ2からＰNに順次変え第２ステップから第２０ステップを繰返して第２等濃度線ＨＣ２Xの変曲点の最小値と第２分布線ＢＣ２Xの変曲点の最大値との差ΔＺ2を求めた後にこの差ΔＺ2が最も大きくなる単結晶製造条件を求める第２１ステップとを含むコンピュータを用いて単結晶の無欠陥領域を最大化するシミュレーション方法である。
【００１３】
この請求項２に記載された単結晶の無欠陥領域を最大化するシミュレーション方法では、パラメータＰ1の単結晶製造条件で上記請求項１と同様に、単結晶１４内のボイドの濃度分布及びサイズ分布を予測するとともに、単結晶１４内の高温酸素析出物の濃度分布及びサイズ分布も予測する。次いで単結晶１４内の第２等濃度線ＨＣ２Xと第２分布線ＢＣ２Xを計算により求めた後に、第２等濃度線ＨＣ２Xの変曲点の最大値と第２分布線ＢＣ２Xの変曲点の最小値との差ΔＺ2を計算により求める。
次に上記単結晶製造条件のパラメータをＰ2から順にＰNまで変え、上記と同様にして第２等濃度線ＨＣ２Xの変曲点の最大値と第２分布線ＢＣ２Xの変曲点の最小値との差ΔＺ2を計算により求める。更に上記第２等濃度線ＨＣ２Xの変曲点の最大値と第２分布線ＢＣ２Xの変曲点の最小値との差ΔＺ2が最も大きくなる単結晶製造条件を求める。これにより単結晶１４の引上げ方向及び半径方向に無欠陥領域が最も拡大する単結晶製造条件を計算により正確に求めることができる。
【００１４】
【発明の実施の形態】
次に本発明の実施の形態を図面に基づいて説明する。
図６に示すように、シリコン単結晶引上げ機１１のチャンバ内には、シリコン融液１２を貯留する石英るつぼ１３が設けられる。この石英るつぼ１３は図示しないが黒鉛サセプタ及び支軸を介してるつぼ駆動手段に接続され、るつぼ駆動手段は石英るつぼ１３を回転させるとともに昇降させるように構成される。また石英るつぼ１３の外周面は石英るつぼ１３から所定の間隔をあけてヒータ（図示せず）により包囲され、このヒータは保温筒（図示せず）により包囲される。ヒータは石英るつぼ１３に投入された高純度のシリコン多結晶体を加熱・溶融してシリコン融液１２にする。またチャンバの上端には図示しないが円筒状のケーシングが接続され、このケーシングには引上げ手段が設けられる。引上げ手段はシリコン単結晶１４を回転させながら引上げるように構成される。
【００１５】
このように構成されたシリコン単結晶引上げ機１１におけるシリコン単結晶１４の無欠陥領域を最大化する第１のシミュレーション方法を図１〜図７に基づいて説明する。
▲１▼ シリコン単結晶１４内のボイドの濃度分布及びサイズ分布
先ず第１ステップとして、シリコン単結晶引上げ機１１によりシリコン単結晶１４を引上げるときの単結晶製造条件を任意に定義する。この単結晶製造条件は引上げ機１１によりシリコン単結晶１１を引上げるときに、後述する引上げ機１１のホットゾーンにフィードバックされる変数を一定の間隔で変量させたパラメータ群Ｐ1，Ｐ2，…，ＰNである。また単結晶製造条件としては、シリコン単結晶１４の引上げ速度、シリコン単結晶１４の回転速度、石英るつぼ１３の回転速度、アルゴンガスの流量、ヒートキャップを構成する部材の形状及び材質、ヒートキャップの下端及びシリコン融液１２表面間のギャップ、ヒータ出力等が挙げられる。
【００１６】
第２ステップとして、シリコン単結晶１４を所定長さＬ1（例えば１００ｍｍ）まで引上げた状態における引上げ機１１のホットゾーンの各部材、即ちチャンバ，石英るつぼ１３，シリコン融液１２，シリコン単結晶１４，黒鉛サセプタ，保温筒等をメッシュ分割してモデル化する。具体的には上記ホットゾーンの各部材のメッシュ点の座標データをコンピュータに入力する。このときシリコン融液１２のメッシュのうちシリコン単結晶１４の径方向のメッシュであってかつシリコン融液１２のシリコン単結晶１４直下の一部又は全部のメッシュ（以下、径方向メッシュという。）を０．０１〜５．００ｍｍ、好ましくは０．２５〜１．００ｍｍに設定する。またシリコン融液１２のメッシュのうちシリコン単結晶１４の長手方向のメッシュであってかつシリコン融液１２の一部又は全部のメッシュ（以下、長手方向メッシュという。）を０．０１〜５．００ｍｍ、好ましくは０．１〜０．５ｍｍに設定する。
【００１７】
径方向メッシュを０．０１〜５．００ｍｍの範囲に限定したのは、０．０１ｍｍ未満では計算時間が極めて長くなり、５．００ｍｍを越えると計算が不安定になり、繰返し計算を行っても固液界面形状が一定に定まらなくなるからである。また長手方向メッシュを０．０１〜５．００ｍｍの範囲に限定したのは、０．０１ｍｍ未満では計算時間が極めて長くなり、５．００ｍｍを越えると固液界面形状の計算値が実測値と一致しなくなるからである。なお、径方向メッシュの一部を０．０１〜５．００の範囲に限定する場合には、シリコン単結晶１４直下のシリコン融液１２のうちシリコン単結晶１４外周縁近傍のシリコン融液１２を上記範囲に限定することが好ましく、長手方向メッシュの一部を０．０１〜５．００の範囲に限定する場合には、シリコン融液１２の液面近傍及び底近傍を上記範囲に限定することが好ましい。
【００１８】
第３ステップとして上記ホットゾーンの各部材毎にメッシュをまとめ、かつこのまとめられたメッシュに対して各部材の物性値をそれぞれコンピュータに入力する。例えば、チャンバがステンレス鋼にて形成されていれば、そのステンレス鋼の熱伝導率，輻射率，粘性率，体積膨張係数，密度及び比熱がコンピュータに入力される。またシリコン単結晶１４の引上げ長及びこの引上げ長に対応するシリコン単結晶１４の引上げ速度と、後述する乱流モデル式（１）の乱流パラメータＣとをコンピュータに入力する。
【００１９】
第４ステップとして、ホットゾーンの各部材の表面温度分布をヒータの発熱量及び各部材の輻射率に基づいてコンピュータを用いて計算により求める。即ち、ヒータの発熱量を任意に設定してコンピュータに入力するとともに、各部材の輻射率から各部材の表面温度分布をコンピュータを用いて計算により求める。次に第５ステップとしてホットゾーンの各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式（１）をコンピュータを用いて解くことにより各部材の内部温度分布を計算により求める。ここでは、記述を簡単にするためｘｙｚ直交座標系を用いたが、実際の計算では円筒座標系を用いる。
【００２０】
【数１】
一方、シリコン融液１２に関しては、上記熱伝導方程式（１）でシリコン融液１２の内部温度分布を求めた後に、このシリコン融液１２の内部温度分布に基づき、シリコン融液１２が乱流であると仮定して得られた乱流モデル式（２）及びナビエ・ストークスの方程式（３）〜（５）を連結して、シリコン融液１２の内部流速分布をコンピュータを用いて求める。
【００２１】
【数２】
【００２２】
【数３】
【００２３】
上記乱流モデル式（２）はｋｌ（ケイエル）−モデル式と呼ばれ、このモデル式の乱流パラメータＣは０．４〜０．６の範囲内の任意の値が用いられることが好ましい。乱流パラメータＣを０．４〜０．６の範囲に限定したのは、０．４未満又は０．６を越えると計算により求めた界面形状が実測値と一致しないという不具合があるからである。また上記ナビエ・ストークスの方程式（３）〜（５）はシリコン融液１２が非圧縮性であって粘度が一定である流体としたときの運動方程式である。
上記求められたシリコン融液１２の内部流速分布に基づいて熱エネルギ方程式（６）を解くことにより、シリコン融液１２の対流を考慮したシリコン融液１２の内部温度分布をコンピュータを用いて更に求める。
【００２４】
【数４】
【００２５】
次いで第６ステップとして、シリコン単結晶１４及びシリコン融液１２の固液界面形状を図６の点Ｓで示すシリコンの三重点Ｓ（固体と液体と気体の三重点(tri-junction)）を含む等温線に合せてコンピュータを用いて計算により求める。第７ステップとして、コンピュータに入力するヒータの発熱量を変更し（次第に増大し）、上記第４ステップから第６ステップを三重点Ｓがシリコン単結晶１４の融点になるまで繰返した後に、引上げ機１１内の温度分布を計算してシリコン単結晶のメッシュの座標及び温度を求め、これらのデータをコンピュータに記憶させる。
【００２６】
次に第８ステップとして、シリコン単結晶１４の引上げ長Ｌ1にδ（例えば５０ｍｍ）だけ加えて上記第２ステップから第７ステップまでを繰返した後に、引上げ機１１内の温度分布を計算してシリコン単結晶１４のメッシュの座標及び温度を求め、これらのデータをコンピュータに記憶させる。この第８ステップは、シリコン単結晶１４の引上げ長Ｌ1が長さＬ2（Ｌ2はシリコン融液１２から切離されたときのシリコン単結晶１４の長さ（成長完了時の結晶長）である。）に達してシリコン単結晶１４がシリコン融液１２から切離された後、更にシリコン単結晶１４が引上げられてその高さＨ1（Ｈ1はシリコン単結晶１４の直胴開始部からシリコン融液１２の液面までの距離である（図６）。）がＨ2（Ｈ2は冷却完了時のシリコン単結晶１４の直胴開始部からシリコン融液１２の液面までの距離である。）に達するまで、即ちシリコン単結晶１４の冷却が完了するまで行われる。なお、シリコン単結晶１４がシリコン融液１２から切離された後は、シリコン単結晶１４の引上げ高さＨ1にδ（例えば５０ｍｍ）だけ加え、上記と同様に上記第２ステップから第７ステップまでを繰返す。
【００２７】
シリコン単結晶１４の引上げ高さＨ1がＨ2に達すると、第９ステップに移行する。第９ステップでは、シリコン単結晶１４をシリコン融液１２から成長させて引上げ始めたときｔ0から、シリコン単結晶１４をシリコン融液１２から切離して更にシリコン単結晶１４を引上げ、その冷却が完了したときｔ1までの時間を、所定の間隔Δｔ秒（微小時間間隔）毎に区切る。このときシリコン単結晶１４内の格子間シリコン及び空孔の拡散係数及び境界条件のみならず、後述するボイド及び高温酸素析出物の濃度分布及びサイズ分布を求めるための式に用いられる定数をそれぞれコンピュータに入力する。上記区切られた時間間隔Δｔ秒毎に、第８ステップで求めたシリコン単結晶１４のメッシュの座標及び温度のデータから、シリコン単結晶１４の引上げ長Ｌ1及び引上げ高さＨ1と、シリコン単結晶１４内の温度分布とを計算により求める。
【００２８】
即ち、第２〜第８ステップでシリコン単結晶のメッシュの座標及び温度を引上げ長δ毎に求め、シリコン単結晶を例えば５０ｍｍ引上げるのに数十分要するため、この数十分間でのシリコン単結晶のメッシュの温度変化を時間の関数として微分することにより、時刻ｔ0からΔｔ秒後におけるシリコン単結晶１４の引上げ長Ｌ1及び引上げ高さＨ1とシリコン単結晶１４内の温度分布を算出する。次にシリコン単結晶１４内の空孔及び格子間シリコンの拡散係数及び境界条件に基づいて拡散方程式を解くことにより、Δｔ秒経過後の空孔及び格子間シリコンの濃度分布を計算により求める（第１０ステップ）。
【００２９】
具体的には、空孔の濃度Ｃvの計算式が次の式（７）で示され、格子間シリコンの濃度Ｃiの計算式が次の式（８）で示される。式（７）及び式（８）において、濃度Ｃv及び濃度Ｃiの経時的進展を計算するために、空孔と格子間シリコンの熱平衡がシリコン単結晶の全表面で維持されると仮定する。
【００３０】
【数５】
上記平衡式は時間で微分され、空孔及び格子間シリコンに対してそれぞれ次の式（９）及び式（１０）になる。
【００３１】
【数６】
【００３２】
【数７】
【数８】
【００３３】
一方、上記式（９）が成立つのは、空孔が析出するための流速が十分大きく、ＳｉマトリックスとＳｉＯ2との単位質量当りの体積差を埋められる場合、即ちＤv（Ｃv−Ｃv e）≧γＤ0Ｃ0の場合である。上記以外の場合には、次式（１３）が成立つ。
【００３４】
【数９】
【００３５】
【数１０】
【００３６】
第１２ステップは後述することにして、第１３ステップを説明する。第１３ステップでは、シリコン単結晶１４内のそれぞれのメッシュの格子点における温度が次第に低下してボイドの形成開始温度Ｔv0になったときに、次の近似式（１５）を用いてボイドの濃度Ｎv0を求める。
【００３７】
【数１１】
【００３８】
【数１２】
【００３９】
▲２▼ シリコン単結晶１４内の高温酸素析出物の濃度分布及びサイズ分布
第１２ステップに戻って、上記拡散方程式（９）及び（１０）を解くことにより求めた空孔の濃度Ｃv分布に基づいて、高温酸素析出物の形成開始温度Ｔpを次の式（１７）から計算により求める。
【００４０】
【数１３】
【００４１】
次に第１５ステップとして、シリコン単結晶１４内のそれぞれのメッシュの格子点における温度が次第に低下して高温酸素析出物の形成開始温度Ｔpになったときに、次の近似式（１８）を用いて高温酸素析出物の濃度Ｎpを求める。なお、式（１８）において、ａ2は定数である。
【００４２】
【数１４】
【００４３】
【数１５】
【００４４】
一方、上記式（１９）が成立つのは、空孔が析出するための流速が十分大きく、ＳｉマトリックスとＳｉＯ2との単位質量当りの体積差を埋められる場合、即ちＤv（Ｃv−Ｃv e）≧γＤ0Ｃ0の場合である。上記以外の場合には、次式（２０）が成立つ。
【００４５】
【数１６】
【００４６】
第１８ステップとして、上記高温酸素析出物の存在するＰバンドとボイド及び高温酸素析出物のいずれも存在しない無欠陥領域とを区画し、かつ無欠陥領域がＰバンドより単結晶ボトム側に位置するように区画する第１等濃度線ＨＣ１Xを計算により求める（図７）。次に第１９ステップとして、高濃度酸素析出物が存在するＢバンドと無欠陥領域とを区画し、かつ無欠陥領域がＢバンドより単結晶トップ側に位置するように区画する第１分布線ＢＣ１Xを計算により求める。上記第１分布線は、格子間シリコン濃度がシリコン融点での格子間シリコンの平衡濃度に対して０．１２〜０．１３％の範囲の特定値、例えば０．１２６％である等濃度線をいう。ここで格子間シリコン濃度０．１２〜０．１３％の範囲の特定値、例えば０．１２６％とは、シミュレーションを行って得られた格子間シリコン濃度のうちライフタイムのマップ等から観察されるＢバンドに対応する濃度である。またＢバンドとは、格子間シリコンの凝集体が核となり熱処理によって酸素析出が高濃度に発生している領域をいう。
【００４７】
その後、第２０ステップとして、第１等濃度線ＨＣ１Xの変曲点の最大値と第１分布線ＢＣ１Xの変曲点の最小値との差ΔＺ1(X)を求める。第１等濃度線ＨＣ１Xの変曲点が例えば３つある場合には、各変曲点Ｑ１(1,X)、Ｑ１(2,X)及びＱ１(3,X)の座標を(ｒQ1(1,X)，ｚQ1(1,X))、(ｒQ1(2,X)，ｚQ1(2,X))及び(ｒQ1(3,X)，ｚQ1(3,X))とすると、これらの座標のうちｚ座標ｚQ1(1,X)、ｚQ1(2,X)及びｚQ1(3,X)のうちの最大値ｚQ1(1,X)を求める。また第１分布線ＢＣ１Xの変曲点が例えば５つある場合には、各変曲点Ｓ１(1,X)、Ｓ１(2,X)、Ｓ１(3,X)、Ｓ１(4,X)及びＳ１(5,X)の座標を(ｒS1(1,X)，ｚS1(1,X))、(ｒS1(2,X)，ｚS1(2,X))、(ｒS1(3,X)，ｚS1(3,X))、(ｒS1(4,X)，ｚS1(4,X))及び(ｒS1(5,X)，ｚS1(5,X))とすると、これらの座標のうちｚ座標ｚS1(1,X)、ｚS1(2,X)、ｚS1(3,X)、ｚS1(4,X)及びｚS1(5,X)のうちの最小値ｚS1(3,X)を求める。そして第１等濃度線ＨＣ１Xの変曲点の最大値ｚQ1(1,X)と第１分布線ＢＣ１Xの変曲点の最小値ｚS1(3,X)との差ΔＺ1(X)を求める。
【００４８】
第２１ステップとして、単結晶製造条件のパラメータをＰ2に変え第２ステップから第２０ステップを実行し、単結晶製造条件のパラメータをＰ3に変え第２ステップから第２０ステップを実行するというように、単結晶製造条件のパラメータがＰNになるまで上記第２ステップから第２０ステップを繰返した後に、第１等濃度線ＨＣ１Xの変曲点の最大値と第１分布線ＢＣ１Xの変曲点の最小値との差ΔＺ1(X)を求め、更にこの差ΔＺ1(X)が最も大きくなる単結晶製造条件を計算により求める。
【００４９】
このように計算して求められたシリコン単結晶引上げ機１１における単結晶製造条件は、シリコン単結晶１４の引上げ方向及び半径方向に無欠陥領域が最も拡大する条件であり、この単結晶製造条件を上記計算により正確に求めることができる。この結果、上述のようにコンピュータを用いた計算によりシリコン単結晶１４の無欠陥領域が最も拡大する単結晶製造条件を設定することができるので、実際にシリコン単結晶を引上げてライフタイム等を測定する実験は上記単結晶製造条件が無欠陥領域が実際に拡大しているか否かの確認のためにだけ行えば済む。従って、上記実験の回数を低減できるので、比較的少ない時間と労力で、最適な単結晶製造条件を得ることができる。なお、図７から明らかなように、上記差ΔＺ1(X)が大きくなると、シリコン単結晶１４の引上げ方向及び半径方向にわたって無欠陥領域とするためのシリコン単結晶１４の引上げ速度の幅、即ちピュアマージンが大きくなるので、全長のわたって無欠陥のシリコン単結晶１４を比較的容易に引上げることができる。
【００５０】
次にシリコン単結晶引上げ機１１におけるシリコン単結晶１４の無欠陥領域を最大化する第２のシミュレーション方法を図６〜図１２に基づいて説明する。
第１〜第１７ステップは、上記第１のシミュレーション方法の第１〜第１７ステップと同一である。
【００５１】
第１８ステップとして、上記高温酸素析出物の存在するＰバンドとボイド及び高温酸素析出物のいずれも存在しない無欠陥領域とを区画し、かつ無欠陥領域がＰバンドより単結晶トップ側に位置するように区画する第２等濃度線ＨＣ２Xを計算により求める（図７）。次に第１９ステップとして、高濃度酸素析出物が存在するＢバンドと無欠陥領域とを区画し、かつ無欠陥領域がＢバンドより単結晶ボトム側に位置するように区画する第２分布線ＢＣ２Xを計算により求める。上記第２分布線は、格子間シリコン濃度がシリコン融点での格子間シリコンの平衡濃度に対して０．１２〜０．１３％の範囲の特定値、例えば０．１２６％である等濃度線をいう。ここで格子間シリコン濃度０．１２〜０．１３％の範囲の特定値、例えば０．１２６％とは、シミュレーションを行って得られた格子間シリコン濃度のうちライフタイムのマップ等から観察されるＢバンドに対応する濃度である。
【００５２】
その後、第２０ステップとして、第２等濃度線ＨＣ２Xの変曲点の最小値と第２分布線ＢＣ２Xの変曲点の最大値との差ΔＺ2(X)を求める。第２等濃度線ＨＣ２Xの変曲点が例えば３つある場合には、各変曲点Ｑ２(1,X)、Ｑ２(2,X)及びＱ２(3,X)の座標を(ｒQ2(1,X)，ｚQ2(1,X))、(ｒQ2(2,X)，ｚQ2(2,X))及び(ｒQ2(3,X)，ｚQ2(3,X))とすると、これらの座標のうちｚ座標ｚQ2(1,X)、ｚQ2(2,X)及びｚQ2(3,X)のうちの最小値ｚQ2(1,X)を求める。また第２分布線ＢＣ２Xの変曲点が例えば５つある場合には、各変曲点Ｓ２(1,X)、Ｓ２(2,X)、Ｓ２(3,X)、Ｓ２(4,X)及びＳ２(5,X)の座標を(ｒS2(1,X)，ｚS2(1,X))、(ｒS2(2,X)，ｚS2(2,X))、(ｒS2(3,X)，ｚS2(3,X))、(ｒS2(4,X)，ｚS2(4,X))及び(ｒS2(5,X)，ｚS2(5,X))とすると、これらの座標のうちｚ座標ｚS2(1,X)、ｚS2(2,X)、ｚS2(3,X)、ｚS2(4,X)及びｚS2(5,X)のうちの最大値ｚS2(3,X)を求める。そして第２等濃度線ＨＣ２Xの変曲点の最小値ｚQ2(1,X)と第２分布線ＢＣ２Xの変曲点の最大値ｚS2(3,X)との差ΔＺ2(X)を求める。
【００５３】
第２１ステップとして、単結晶製造条件のパラメータをＰ2に変え第２ステップから第２０ステップを実行し、単結晶製造条件のパラメータをＰ3に変え第２ステップから第２０ステップを実行するというように、単結晶製造条件のパラメータがＰNになるまで上記第２ステップから第２０ステップを繰返した後に、第２等濃度線ＨＣ２Xの変曲点の最小値と第２分布線ＢＣ２Xの変曲点の最大値との差ΔＺ2(X)を求め、更にこの差ΔＺ2(X)が最も大きくなる単結晶製造条件を計算により求める。
【００５４】
このように計算して求められたシリコン単結晶引上げ機１１における単結晶製造条件は、シリコン単結晶１４の引上げ方向及び半径方向に無欠陥領域が最も拡大する条件であり、この単結晶製造条件を上記計算により正確に求めることができる。この結果、上述のようにコンピュータを用いた計算によりシリコン単結晶１４の無欠陥領域が最も拡大する単結晶製造条件を設定することができるので、実際にシリコン単結晶を引上げてライフタイム等を測定する実験は上記単結晶製造条件が無欠陥領域が実際に拡大しているか否かの確認のためにだけ行えば済む。従って、上記実験の回数を低減できるので、比較的少ない時間と労力で、最適な単結晶製造条件を得ることができる。なお、図７から明らかなように、上記差ΔＺ2(X)が大きくなると、シリコン単結晶１４の引上げ方向及び半径方向にわたって無欠陥領域とするためのシリコン単結晶１４の引上げ速度の幅、即ちピュアマージンが大きくなるので、全長のわたって無欠陥のシリコン単結晶１４を比較的容易に引上げることができる。
なお、上記実施の形態では、単結晶としてシリコン単結晶を挙げたが、ＧａＡｓ単結晶，ＩｎＰ単結晶，ＺｎＳ単結晶若しくはＺｎＳｅ単結晶でもよい。
【００５５】
【発明の効果】
以上述べたように、本発明によれば、先ず所定の単結晶製造条件で融液の対流を考慮することにより、計算した結晶及び融液の固液界面形状が実際の単結晶の引上げ時の形状とほぼ一致するようにして単結晶の内部温度を求めるだけでなく、更に冷却過程における単結晶内の温度分布までも求めることによって、即ち融液から切離された単結晶の冷却過程における単結晶の徐冷及び急冷の効果を考慮することによって、単結晶内の欠陥の濃度分布及びサイズ分布をコンピュータを用いて求める。次に単結晶内の第１等濃度線と第１分布線を計算により求めて、第１等濃度線の変曲点の最大値と第１分布線の変曲点の最小値との差を計算により求めた後に、上記単結晶製造条件のパラメータを変え、上記と同様にして第１等濃度線の変曲点の最大値と第１分布線の変曲点の最小値との差を計算により求める。更に上記第１等濃度線の変曲点の最大値と第１分布線の変曲点の最小値との差が最も大きくなる単結晶製造条件を求める。これにより単結晶の引上げ方向及び半径方向に無欠陥領域が最も拡大する単結晶製造条件を計算により正確に求めることができる。このため上記無欠陥領域が最も拡大する単結晶製造条件を設定するのに数多くの実験を繰返す必要のあった従来の無欠陥の単結晶インゴットと比較して、本発明では、上記無欠陥領域が最も拡大する単結晶製造条件を設定するための実験の回数を低減できる。この結果、本発明では、比較的少ない時間と労力で、無欠陥領域が最も拡大する単結晶製造条件を得ることができる。
【００５６】
また単結晶内の第２等濃度線と第２分布線を計算により求めて、第２等濃度線の変曲点の最大値と第２分布線の変曲点の最小値との差を計算により求めた後に、上記単結晶製造条件のパラメータを変え、上記と同様にして第２等濃度線の変曲点の最大値と第２分布線の変曲点の最小値との差を計算により求め、更に上記第１等濃度線の変曲点の最大値と第１分布線の変曲点の最小値との差が最も大きくなる単結晶製造条件を求めても、上記と同様の効果が得られる。
【図面の簡単な説明】
【図１】本発明実施形態の第１のシミュレーション方法の第１段を示すフローチャート。
【図２】そのシミュレーション方法の第２段を示すフローチャート。
【図３】そのシミュレーション方法の第３段を示すフローチャート。
【図４】そのシミュレーション方法の第４段を示すフローチャート。
【図５】そのシミュレーション方法の第５段を示すフローチャート。
【図６】シリコン融液をメッシュ構造としたシリコン単結晶の引上げ機の要部断面図。
【図７】そのシリコン単結晶の引上げ速度を変化させたときのシリコン単結晶内の格子間シリコン及び空孔の分布を示す説明図。
【図８】本発明実施形態の第２のシミュレーション方法の第１段を示すフローチャート。
【図９】そのシミュレーション方法の第２段を示すフローチャート。
【図１０】そのシミュレーション方法の第３段を示すフローチャート。
【図１１】そのシミュレーション方法の第４段を示すフローチャート。
【図１２】そのシミュレーション方法の第５段を示すフローチャート。
【符号の説明】
１１ シリコン単結晶引上げ機
１２ シリコン融液
１４ シリコン単結晶
Ｓ シリコンの三重点[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a method for maximizing a defect-free region of a single crystal such as silicon pulled by a Czochralski (hereinafter referred to as CZ) method by computer simulation.
[0002]
[Prior art]
Factors that reduce device yields due to the recent miniaturization of semiconductor integrated circuits include crystal-origin particles (hereinafter referred to as COP), oxidation-induced stacking faults (hereinafter referred to as COP), The presence of microdefects in oxygen precipitates, which are the core of OISF, or interstitial-type large dislocation (hereinafter referred to as L / D).
[0003]
COP is a pit that appears on the wafer surface when a mirror-polished silicon wafer is washed with a mixed solution of ammonia and hydrogen peroxide. When this wafer is measured with a particle counter, the pits are detected as particles. This pit is caused by crystals. COP causes deterioration of electrical characteristics, for example, dielectric breakdown characteristics (Time Dependent dielectric Breakdown, TDDB) of oxide film, breakdown voltage characteristics of oxide film (Time Zero Dielectric Breakdown, TZDB), and the like. Further, if COP exists on the wafer surface, a step is generated in the device wiring process, which may cause disconnection. In addition, the element isolation portion also causes leakage and the like, thereby reducing the product yield.
[0004]
OISF is a defect that is manifested in a thermal oxidation process or the like when a semiconductor device is manufactured by introducing minute oxygen precipitation nuclei into a silicon single crystal during crystal growth. This OISF causes a failure such as increasing the leakage current of the device. L / D is also referred to as a dislocation cluster, or it is also referred to as a dislocation pit because a pit is generated when a silicon wafer having such a defect is immersed in a selective etching solution containing hydrofluoric acid as a main component. This L / D also causes deterioration of electrical characteristics such as leakage characteristics and isolation characteristics. From the above, it is necessary to reduce COP, OISF, and L / D from a silicon wafer used for manufacturing a semiconductor integrated circuit.
[0005]
A defect-free single crystal ingot having no COP, OISF, and L / D and a silicon wafer sliced from the single crystal ingot are disclosed (for example, see Patent Document 1). This defect-free single crystal ingot has a perfect region [P] when a perfect region in which no agglomerates of vacancy-type point defects and agglomerates of interstitial silicon type point defects exist in the single crystal ingot, respectively. P] is a single crystal ingot. The perfect region [P] includes a region [V] having point defects in which vacancy-type point defects are dominant and agglomerated in a single crystal ingot, and a point defect in which interstitial silicon type point defects are dominant and agglomerated. It is interposed between the region [I].
[0006]
[Patent Document 1]
Japanese Patent Laid-Open No. 11-1393 corresponding to US Pat. No. 6,045,610
[0007]
[Problems to be solved by the invention]
However, the conventional defect-free single crystal ingot disclosed in Japanese Patent Application Laid-Open No. 11-1393 has a void-type point defect aggregate and an interstitial silicon-type point defect aggregate over the entire length of the single crystal ingot. In order to set the pulling conditions that will not be performed, it is necessary to repeat a number of experiments in which the single crystal ingot is actually pulled and the lifetime measured for each type of puller and the thickness of the single crystal ingot. There was a problem that required a great deal of time and effort to set up.
In addition, in the conventional defect-free single crystal ingot, the V / G theory used for obtaining the point defect concentration, that is, the pulling rate of the single crystal is V (mm / min), and the solid crystals and the melt are solidified. When the temperature gradient near the liquid interface is G (° C./mm), V / G (mm2Because it is a one-dimensional model (model only in the pulling direction of the single crystal) that simplifies the Boronkov theory that controls (/ min · ° C), the point defect diffusion (outward diffusion) in the radial direction of the single crystal is considered. However, there is a problem that the concentration gradient of point defects in the radial direction of the single crystal cannot be obtained.
Further, in the conventional defect-free single crystal ingot, when the concentration of point defects is determined by V / G theory, it is necessary to model and calculate the convection of the melt in order to determine G with high accuracy. Nevertheless, since such consideration is not made at all, there is a problem that the concentration distribution of point defects in the single crystal cannot be obtained with high accuracy.
[0008]
The object of the present invention is to reduce the number of experiments of actually pulling up a single crystal, and thereby the defect-free region is expanded most in the pulling direction and radial direction of the single crystal with relatively little time and effort. An object of the present invention is to provide a simulation method capable of accurately obtaining single crystal production conditions by calculation.
[0009]
[Means for Solving the Problems]
As shown in FIGS. 1 to 7, the invention according to
[0010]
In the simulation method for maximizing the defect-free region of the single crystal according to
Next, the parameters of the above single crystal manufacturing conditions are set to P2P in order fromNUntil the first isoconcentration line HC1XMaximum inflection point and first distribution line BC1XThe difference ΔZ from the minimum value of the inflection point1Is calculated. Further, the first isoconcentration line HC1XMaximum inflection point and first distribution line BC1XThe difference ΔZ from the minimum value of the inflection point1The condition for producing a single crystal with the largest value is obtained. Thereby, the single crystal manufacturing conditions under which the defect-free region is expanded most in the pulling direction and the radial direction of the
[0011]
In the present specification, the “high temperature oxygen precipitate” refers to a material in which the inner wall oxide film precipitated inside the void defect is treated like a sphere. The reason why the oxygen precipitate is limited to a high temperature is to distinguish it from a precipitate nucleus such as BMD (Bulk Micro Defect) that precipitates at a temperature lower than the void defect. In this specification, the “defect-free region” includes a defect-free region including an intrinsic point defect in which vacancies are dominant and a defect-free region including an intrinsic defect in which interstitial atoms are dominant.
[0012]
As shown in FIG. 6 and FIGS. 8 to 12, the invention according to
[0013]
In the simulation method for maximizing the defect-free region of the single crystal according to
Next, the parameters of the above single crystal manufacturing conditions are set to P2P in order fromNUntil the second isoconcentration line HC2XMaximum inflection point and second distribution line BC2XThe difference ΔZ from the minimum value of the inflection point2Is calculated. Further, the second isoconcentration line HC2XMaximum inflection point and second distribution line BC2XThe difference ΔZ from the minimum value of the inflection point2The condition for producing a single crystal with the largest value is obtained. Thereby, the single crystal manufacturing conditions under which the defect-free region is expanded most in the pulling direction and the radial direction of the
[0014]
DETAILED DESCRIPTION OF THE INVENTION
Next, embodiments of the present invention will be described with reference to the drawings.
As shown in FIG. 6, a
[0015]
A first simulation method for maximizing the defect-free region of the silicon
(1) Concentration distribution and size distribution of voids in the silicon
First, as a first step, single crystal manufacturing conditions for pulling the silicon
[0016]
As a second step, the silicon
[0017]
The reason why the radial mesh is limited to the range of 0.01 to 5.00 mm is that the calculation time becomes extremely long if it is less than 0.01 mm, and the calculation becomes unstable if it exceeds 5.00 mm. This is because the solid-liquid interface shape cannot be fixed. The longitudinal mesh is limited to the range of 0.01 to 5.00 mm because the calculation time is extremely long if it is less than 0.01 mm, and the calculated value of the solid-liquid interface shape matches the actual measurement value if it exceeds 5.00 mm. Because it will not do. When a part of the radial mesh is limited to the range of 0.01 to 5.00, the
[0018]
As a third step, the meshes are grouped for each member in the hot zone, and the physical property values of the members are input to the computer for the grouped meshes. For example, if the chamber is made of stainless steel, the thermal conductivity, emissivity, viscosity, volume expansion coefficient, density and specific heat of the stainless steel are input to the computer. Further, the pulling length of the silicon
[0019]
As a fourth step, the surface temperature distribution of each member in the hot zone is obtained by calculation using a computer based on the amount of heat generated by the heater and the radiation rate of each member. That is, the heating value of the heater is arbitrarily set and input to the computer, and the surface temperature distribution of each member is obtained by calculation using the computer from the emissivity of each member. Next, as a fifth step, the internal temperature distribution of each member is calculated by solving the heat conduction equation (1) using a computer based on the surface temperature distribution and the thermal conductivity of each member in the hot zone. Here, in order to simplify the description, an xyz orthogonal coordinate system is used, but a cylindrical coordinate system is used in actual calculation.
[0020]
[Expression 1]
On the other hand, for the
[0021]
[Expression 2]
[0022]
[Equation 3]
[0023]
The turbulent model equation (2) is referred to as a kl-model equation, and an arbitrary value within the range of 0.4 to 0.6 is preferably used as the turbulent parameter C of the model equation. The reason why the turbulent flow parameter C is limited to the range of 0.4 to 0.6 is that when the value is less than 0.4 or exceeds 0.6, there is a problem that the interface shape obtained by calculation does not match the actual measurement value. . The Navier-Stokes equations (3) to (5) are equations of motion when the
By solving the thermal energy equation (6) based on the obtained internal flow velocity distribution of the
[0024]
[Expression 4]
[0025]
Next, as a sixth step, the silicon-
[0026]
Next, as an eighth step, the pulling length L of the silicon
[0027]
Lifting height H of silicon
[0028]
That is, in the second to eighth steps, the coordinates and temperature of the silicon single crystal mesh are obtained for each pulling length δ, and it takes several tens of minutes to pull the silicon single crystal, for example, 50 mm. By differentiating the temperature change of the single crystal mesh as a function of time, the time t0The pulling length L of the silicon
[0029]
Specifically, the vacancy concentration CvIs calculated by the following equation (7), and the interstitial silicon concentration CiIs calculated by the following equation (8). In formula (7) and formula (8), concentration CvAnd concentration CiIn order to calculate the evolution over time, it is assumed that the thermal equilibrium between vacancies and interstitial silicon is maintained on the entire surface of the silicon single crystal.
[0030]
[Equation 5]
The above equilibrium equation is differentiated with respect to time and becomes the following equations (9) and (10) for the vacancies and interstitial silicon, respectively.
[0031]
[Formula 6]
[0032]
[Expression 7]
[Equation 8]
[0033]
On the other hand, the above formula (9) is satisfied because the flow velocity for the precipitation of the vacancies is sufficiently large, and the Si matrix and SiO2And the volume difference per unit mass can be filled, that is, Dv(Cv-Cv e) ≧ γD0C0This is the case. In cases other than the above, the following equation (13) is established.
[0034]
[Equation 9]
[0035]
[Expression 10]
[0036]
The twelfth step will be described later, and the thirteenth step will be described. In the thirteenth step, the temperature at the lattice point of each mesh in the silicon
[0037]
[Expression 11]
[0038]
[Expression 12]
[0039]
(2) Concentration distribution and size distribution of high-temperature oxygen precipitates in the silicon
Returning to the twelfth step, the vacancy concentration C determined by solving the diffusion equations (9) and (10).vBased on the distribution, the formation start temperature T of the high-temperature oxygen precipitatepIs calculated from the following equation (17).
[0040]
[Formula 13]
[0041]
Next, as a fifteenth step, the temperature at the lattice point of each mesh in the silicon
[0042]
[Expression 14]
[0043]
[Expression 15]
[0044]
On the other hand, the above equation (19) is satisfied because the flow rate for the precipitation of vacancies is sufficiently large.2And the volume difference per unit mass can be filled, that is, Dv(Cv-Cv e) ≧ γD0C0This is the case. In cases other than the above, the following equation (20) is established.
[0045]
[Expression 16]
[0046]
As an eighteenth step, the P band where the high-temperature oxygen precipitates are present and a defect-free region where neither voids nor high-temperature oxygen precipitates exist are defined, and the defect-free region is located on the single crystal bottom side from the P band. First concentration line HC1XIs obtained by calculation (FIG. 7). Next, as a nineteenth step, a first distribution line BC1 that partitions a B band in which high-concentration oxygen precipitates are present from a defect-free region and that the defect-free region is positioned on the single crystal top side from the B band.XIs calculated. The first distribution line is an isoconcentration line in which the interstitial silicon concentration is a specific value in the range of 0.12 to 0.13% relative to the equilibrium concentration of interstitial silicon at the silicon melting point, for example, 0.126%. Say. Here, the specific value in the range of the interstitial silicon concentration of 0.12 to 0.13%, for example, 0.126%, is observed from the lifetime map or the like among the interstitial silicon concentrations obtained by performing the simulation. The density corresponds to the B band. The B band refers to a region in which an aggregate of interstitial silicon serves as a nucleus and oxygen precipitation occurs at a high concentration by heat treatment.
[0047]
Thereafter, as the 20th step, the first isoconcentration line HC1XMaximum inflection point and first distribution line BC1XThe difference ΔZ from the minimum value of the inflection point1Find (X). First isoconcentration line HC1XIf there are, for example, three inflection points, the coordinates of each inflection point Q1 (1, X), Q1 (2, X) and Q1 (3, X) are expressed as (rQ1(1, X), zQ1(1, X)), (rQ1(2, X), zQ1(2, X)) and (rQ1(3, X), zQ1(3, X)) of these coordinates, z coordinate zQ1(1, X), zQ1(2, X) and zQ1The maximum value z of (3, X)Q1Find (1, X). The first distribution line BC1XFor example, when there are five inflection points, each inflection point S1 (1, X), S1 (2, X), S1 (3, X), S1 (4, X) and S1 (5, X) ) Coordinate (rS1(1, X), zS1(1, X)), (rS1(2, X), zS1(2, X)), (rS1(3, X), zS1(3, X)), (rS1(4, X), zS1(4, X)) and (rS1(5, X), zS1(5, X)) of these coordinates, z coordinate zS1(1, X), zS1(2, X), zS1(3, X), zS1(4, X) and zS1Minimum value z of (5, X)S1Find (3, X). And the first isoconcentration line HC1XMaximum value of inflection point zQ1(1, X) and first distribution line BC1XMinimum value z of inflection pointS1ΔZ difference from (3, X)1Find (X).
[0048]
As the 21st step, the parameters of the single crystal manufacturing conditions are set to P2Execute
[0049]
The single crystal manufacturing conditions in the silicon single
[0050]
Next, a second simulation method for maximizing the defect-free region of the silicon
The first to seventeenth steps are the same as the first to seventeenth steps of the first simulation method.
[0051]
As an eighteenth step, the P band where the high-temperature oxygen precipitates are present and the defect-free region where neither voids nor high-temperature oxygen precipitates exist are defined, and the defect-free region is located on the single crystal top side from the P band. Second isoconcentration line HC2XIs obtained by calculation (FIG. 7). Next, as a nineteenth step, a second distribution line BC2 that partitions the B band in which high-concentration oxygen precipitates are present and the defect-free region, and partitions the defect-free region so that it is located on the single crystal bottom side from the B band.XIs calculated. The second distribution line is an isoconcentration line in which the interstitial silicon concentration is a specific value in the range of 0.12 to 0.13% relative to the equilibrium concentration of interstitial silicon at the silicon melting point, for example, 0.126%. Say. Here, the specific value in the range of the interstitial silicon concentration of 0.12 to 0.13%, for example, 0.126%, is observed from the lifetime map or the like among the interstitial silicon concentrations obtained by performing the simulation. The density corresponds to the B band.
[0052]
Thereafter, as the 20th step, the second isoconcentration line HC2XMinimum inflection point and second distribution line BC2XThe difference ΔZ from the maximum inflection point2Find (X). Second isoconcentration line HC2XIf there are, for example, three inflection points, the coordinates of each inflection point Q2 (1, X), Q2 (2, X) and Q2 (3, X) are expressed as (rQ2(1, X), zQ2(1, X)), (rQ2(2, X), zQ2(2, X)) and (rQ2(3, X), zQ2(3, X)) of these coordinates, z coordinate zQ2(1, X), zQ2(2, X) and zQ2Minimum value z of (3, X)Q2Find (1, X). The second distribution line BC2XFor example, when there are five inflection points, each inflection point S2 (1, X), S2 (2, X), S2 (3, X), S2 (4, X) and S2 (5, X ) Coordinate (rS2(1, X), zS2(1, X)), (rS2(2, X), zS2(2, X)), (rS2(3, X), zS2(3, X)), (rS2(4, X), zS2(4, X)) and (rS2(5, X), zS2(5, X)) of these coordinates, z coordinate zS2(1, X), zS2(2, X), zS2(3, X), zS2(4, X) and zS2The maximum value z of (5, X)S2Find (3, X). And the second isoconcentration line HC2XMinimum value z of inflection pointQ2(1, X) and second distribution line BC2XMaximum value of inflection point zS2ΔZ difference from (3, X)2Find (X).
[0053]
As the 21st step, the parameters of the single crystal manufacturing conditions are set to P2Execute
[0054]
The single crystal manufacturing conditions in the silicon single
In the above embodiment, a silicon single crystal is used as the single crystal. However, a GaAs single crystal, an InP single crystal, a ZnS single crystal, or a ZnSe single crystal may be used.
[0055]
【The invention's effect】
As described above, according to the present invention, by first considering the convection of the melt under predetermined single crystal manufacturing conditions, the calculated crystal and the solid-liquid interface shape of the melt can be obtained when the actual single crystal is pulled up. In addition to determining the internal temperature of the single crystal so that it substantially matches the shape, it is also possible to determine the temperature distribution within the single crystal during the cooling process, that is, the single crystal in the cooling process of the single crystal separated from the melt. By considering the effect of slow cooling and rapid cooling of the crystal, the concentration distribution and size distribution of defects in the single crystal are obtained using a computer. Next, the first isoconcentration line and the first distribution line in the single crystal are obtained by calculation, and the difference between the maximum value of the inflection point of the first isoconcentration line and the minimum value of the inflection point of the first distribution line is calculated. After obtaining by calculation, the parameters of the single crystal production conditions are changed, and the difference between the maximum value of the inflection point of the first isoconcentration line and the minimum value of the inflection point of the first distribution line is calculated in the same manner as described above. Ask for. Further, a single crystal production condition is obtained in which the difference between the maximum value of the inflection point of the first isoconcentration line and the minimum value of the inflection point of the first distribution line is the largest. As a result, the single crystal manufacturing conditions under which the defect-free region expands most in the pulling direction and the radial direction of the single crystal can be accurately obtained by calculation. For this reason, compared with a conventional defect-free single crystal ingot that required many experiments to set a single crystal production condition in which the defect-free region is most expanded, in the present invention, the defect-free region is The number of experiments for setting the most expanded single crystal production conditions can be reduced. As a result, in the present invention, it is possible to obtain a single crystal manufacturing condition in which the defect-free region is most expanded with relatively little time and effort.
[0056]
Also, the second isoconcentration line and the second distribution line in the single crystal are obtained by calculation, and the difference between the maximum value of the inflection point of the second isoconcentration line and the minimum value of the inflection point of the second distribution line is calculated. Then, the parameter of the single crystal production condition is changed, and the difference between the maximum value of the inflection point of the second isoconcentration line and the minimum value of the inflection point of the second distribution line is calculated in the same manner as above. Further, the same effect as described above can be obtained by obtaining a single crystal production condition in which the difference between the maximum value of the inflection point of the first isoconcentration line and the minimum value of the inflection point of the first distribution line is maximized. can get.
[Brief description of the drawings]
FIG. 1 is a flowchart showing a first stage of a first simulation method according to an embodiment of the present invention.
FIG. 2 is a flowchart showing a second stage of the simulation method.
FIG. 3 is a flowchart showing a third stage of the simulation method.
FIG. 4 is a flowchart showing a fourth stage of the simulation method.
FIG. 5 is a flowchart showing a fifth stage of the simulation method.
FIG. 6 is a cross-sectional view of an essential part of a silicon single crystal pulling machine having a silicon melt mesh structure.
FIG. 7 is an explanatory diagram showing the distribution of interstitial silicon and vacancies in the silicon single crystal when the pulling rate of the silicon single crystal is changed.
FIG. 8 is a flowchart showing a first stage of a second simulation method according to the embodiment of the present invention.
FIG. 9 is a flowchart showing a second stage of the simulation method.
FIG. 10 is a flowchart showing a third stage of the simulation method.
FIG. 11 is a flowchart showing a fourth stage of the simulation method.
FIG. 12 is a flowchart showing the fifth stage of the simulation method.
[Explanation of symbols]
11 Silicon single crystal pulling machine
12 Silicon melt
14 Silicon single crystal
S Triple point of silicon
Claims (2)
前記パラメータ群Ｐ1，Ｐ2，…，ＰNのうちパラメータＰ1を前記単結晶製造条件として前記単結晶(14)の融液からの引上げ開始時から前記単結晶(14)の冷却完了時までの前記引上げ機(11)のホットゾーンをメッシュ構造でモデル化する第２ステップと、
前記ホットゾーンの各部材毎にメッシュをまとめかつこのまとめられたメッシュに対する前記各部材の物性値とともに前記パラメータ群Ｐ1，Ｐ2，…，ＰNのうちのパラメータＰ1をコンピュータに入力する第３ステップと、
前記各部材の表面温度分布をヒータの発熱量及び前記各部材の輻射率に基づいて求める第４ステップと、
前記各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式を解くことにより前記各部材の内部温度分布を求めた後に融液(12)が乱流であると仮定して得られた乱流モデル式及びナビエ・ストークスの方程式を連結して解くことにより対流を考慮した前記融液(12)の内部温度分布を更に求める第５ステップと、
前記単結晶(14)及び前記融液(12)の固液界面形状を前記単結晶の三重点(S)を含む等温線に合せて求める第６ステップと、
前記第４ステップから前記第６ステップを前記三重点(S)が前記単結晶(14)の融点になるまで繰返し前記引上げ機(11)内の温度分布を計算して前記単結晶(14)のメッシュの座標及び温度を求めこれらのデータをそれぞれ前記コンピュータに入力する第７ステップと、
前記単結晶(14)の引上げ長及び引上げ高さを段階的に変えて前記第２ステップから前記第７ステップまでを繰返し前記引上げ機(11)内の温度分布を計算して前記単結晶(14)のメッシュの座標及び温度を求めこれらのデータをそれぞれ前記コンピュータに入力する第８ステップと、
前記単結晶(14)の前記融液(12)からの引上げ開始時から前記単結晶(14)の冷却完了時までの時間を所定の間隔毎に区切り前記区切られた時間間隔毎に第８ステップで求めた前記単結晶(14)のメッシュの座標及び温度のデータから前記単結晶(14)の引上げ長及び引上げ高さと前記単結晶(14)内の温度分布とを求める第９ステップと、
前記単結晶(14)内の空孔及び格子間原子の拡散係数及び境界条件に基づいて拡散方程式を解くことにより前記所定の時間間隔の経過した後の空孔及び格子間原子の濃度分布を求める第１０ステップと、
前記空孔の濃度分布に基づいてボイドの形成開始温度を求める第１１ステップと、
前記空孔の濃度分布に基づいて高温酸素析出物の形成開始温度を求める第１２ステップと、
前記単結晶(14)内のそれぞれのメッシュの格子点における温度が次第に低下して前記ボイドの形成開始温度になったときの前記ボイドの濃度を求める第１３ステップと、
前記単結晶(14)内のそれぞれのメッシュの格子点における温度が前記ボイドの形成開始温度より低いときの前記ボイドの半径を求める第１４ステップと、
前記単結晶(14)内のそれぞれのメッシュの格子点における温度が次第に低下して前記高温酸素析出物の形成開始温度になったときの前記高温酸素析出物の濃度を求める第１５ステップと、
前記単結晶(14)内のそれぞれのメッシュの格子点における温度が前記高温酸素析出物の形成開始温度より低いときの前記高温酸素析出物の半径を求める第１６ステップと、
第９ステップから第１６ステップを前記単結晶(14)の冷却が完了するまで繰返す第１７ステップと、
前記高温酸素析出物の存在するＰバンドと前記ボイド及び前記高温酸素析出物のいずれも存在しない無欠陥領域とを区画しかつ前記無欠陥領域が前記Ｐバンドより単結晶ボトム側に位置するように区画する第１等濃度線(HC1X)を求める第１８ステップと、
高濃度酸素析出物が存在するＢバンドと前記無欠陥領域とを区画しかつ前記無欠陥領域が前記Ｂバンドより単結晶トップ側に位置するように区画する第１分布線(BC1X)を求める第１９ステップと、
前記第１等濃度線(HC1X)の変曲点の最大値と前記第１分布線(BC1X)の変曲点の最小値との差(ΔZ1)を求める第２０ステップと、
前記単結晶製造条件のパラメータをＰ2からＰNに順次変え前記第２ステップから前記第２０ステップを繰返して前記第１等濃度線(HC1X)の変曲点の最大値と前記第１分布線(BC1X)の変曲点の最小値との差(ΔZ1)を求めた後に前記差(ΔZ1)が最も大きくなる単結晶製造条件を求める第２１ステップと
を含むコンピュータを用いて単結晶の無欠陥領域を最大化するシミュレーション方法。Parameter groups P 1 , P 2 ,..., P N are arbitrarily defined by varying variables fed back as a single crystal production condition when pulling up the single crystal 14 by the puller 11. The first step;
From the parameter group P 1 , P 2 ,..., P N , the parameter P 1 is used as the single crystal production condition, and the cooling of the single crystal (14) is completed from the start of the pulling of the single crystal (14) from the melt. A second step of modeling the hot zone of the puller (11) up to a mesh structure;
The summarized mesh for each member of the hot zone and the with physical properties of each member for the gathered mesh parameter groups P 1, P 2, ..., first enter the parameters P 1 of the P N to the computer 3 steps,
A fourth step of obtaining the surface temperature distribution of each member based on the heat value of the heater and the radiation rate of each member;
The turbulence obtained by assuming that the melt (12) is turbulent after obtaining the internal temperature distribution of each member by solving the heat conduction equation based on the surface temperature distribution and thermal conductivity of each member. A fifth step of further determining the internal temperature distribution of the melt (12) considering convection by connecting and solving the flow model equation and the Navier-Stokes equations;
A sixth step of obtaining a solid-liquid interface shape of the single crystal (14) and the melt (12) according to an isotherm including a triple point (S) of the single crystal;
The fourth step to the sixth step are repeated until the triple point (S) reaches the melting point of the single crystal (14), and the temperature distribution in the puller (11) is calculated to calculate the single crystal (14). A seventh step of obtaining the coordinates and temperature of the mesh and inputting these data to the computer, respectively;
The pulling length and the pulling height of the single crystal (14) are changed stepwise, and the temperature distribution in the pulling machine (11) is calculated repeatedly from the second step to the seventh step to calculate the single crystal (14 8) to obtain the coordinates and temperature of the mesh and input these data to the computer,
Eighth step of dividing the time from the start of pulling of the single crystal (14) from the melt (12) to the completion of cooling of the single crystal (14) at predetermined intervals. A ninth step of determining the pulling length and height of the single crystal (14) and the temperature distribution in the single crystal (14) from the mesh coordinates and temperature data of the single crystal (14) determined in
The concentration distribution of vacancies and interstitial atoms after the predetermined time interval has been obtained by solving the diffusion equation based on the diffusion coefficients and boundary conditions of vacancies and interstitial atoms in the single crystal (14) 10th step;
An eleventh step of determining a void formation start temperature based on the concentration distribution of the pores;
A twelfth step of determining a formation start temperature of the high-temperature oxygen precipitate based on the concentration distribution of the pores;
A thirteenth step of determining the concentration of the void when the temperature at the lattice point of each mesh in the single crystal (14) gradually decreases to the formation start temperature of the void;
A fourteenth step of determining a radius of the void when a temperature at a lattice point of each mesh in the single crystal (14) is lower than a formation start temperature of the void;
A fifteenth step of determining the concentration of the high-temperature oxygen precipitates when the temperature at the lattice points of each mesh in the single crystal (14) gradually decreases to reach the formation start temperature of the high-temperature oxygen precipitates;
A sixteenth step of obtaining a radius of the high-temperature oxygen precipitate when the temperature at the lattice point of each mesh in the single crystal (14) is lower than the formation start temperature of the high-temperature oxygen precipitate;
A seventeenth step of repeating the ninth to sixteenth steps until the cooling of the single crystal (14) is completed;
The P band where the high-temperature oxygen precipitates are present and the defect-free region where neither the void nor the high-temperature oxygen precipitates exist are partitioned, and the defect-free region is located on the single crystal bottom side from the P band. An eighteenth step for determining a first isoconcentration line (HC1 X ) to be partitioned;
A first distribution line (BC1 X ) is obtained that partitions the B band in which high-concentration oxygen precipitates are present and the defect-free region, and partitions the defect-free region so that it is located on the single crystal top side from the B band. 19th step;
A twentieth step of obtaining a difference (ΔZ 1 ) between the maximum value of the inflection point of the first isoconcentration line (HC1 X ) and the minimum value of the inflection point of the first distribution line (BC1 X );
The parameter of the single crystal production condition is sequentially changed from P 2 to PN, and the second step to the twentieth step are repeated, so that the maximum value of the inflection point of the first isoconcentration line (HC1 X ) and the first distribution are obtained. A computer using a computer including a twenty-first step of obtaining a single crystal production condition in which the difference (ΔZ 1 ) is largest after obtaining a difference (ΔZ 1 ) from the minimum value of the inflection point of the line (BC1 X ). A simulation method that maximizes the defect-free region of a crystal.
前記パラメータ群Ｐ1，Ｐ2，…，ＰNのうちパラメータＰ1を前記単結晶製造条件として前記単結晶(14)の融液からの引上げ開始時から前記単結晶(14)の冷却完了時までの前記引上げ機(11)のホットゾーンをメッシュ構造でモデル化する第２ステップと、
前記ホットゾーンの各部材毎にメッシュをまとめかつこのまとめられたメッシュに対する前記各部材の物性値とともに前記パラメータ群Ｐ1，Ｐ2，…，ＰNのうちのパラメータＰ1をコンピュータに入力する第３ステップと、
前記各部材の表面温度分布をヒータの発熱量及び前記各部材の輻射率に基づいて求める第４ステップと、
前記各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式を解くことにより前記各部材の内部温度分布を求めた後に融液(12)が乱流であると仮定して得られた乱流モデル式及びナビエ・ストークスの方程式を連結して解くことにより対流を考慮した前記融液(12)の内部温度分布を更に求める第５ステップと、
前記単結晶(14)及び前記融液(12)の固液界面形状を前記単結晶の三重点(S)を含む等温線に合せて求める第６ステップと、
前記第４ステップから前記第６ステップを前記三重点(S)が前記単結晶(14)の融点になるまで繰返し前記引上げ機(11)内の温度分布を計算して前記単結晶(14)のメッシュの座標及び温度を求めこれらのデータをそれぞれ前記コンピュータに入力する第７ステップと、
前記単結晶(14)の引上げ長及び引上げ高さを段階的に変えて前記第２ステップから前記第７ステップまでを繰返し前記引上げ機(11)内の温度分布を計算して前記単結晶(14)のメッシュの座標及び温度を求めこれらのデータをそれぞれ前記コンピュータに入力する第８ステップと、
前記単結晶(14)の前記融液(12)からの引上げ開始時から前記単結晶(14)の冷却完了時までの時間を所定の間隔毎に区切り前記区切られた時間間隔毎に第８ステップで求めた前記単結晶(14)のメッシュの座標及び温度のデータから前記単結晶(14)の引上げ長及び引上げ高さと前記単結晶(14)内の温度分布とを求める第９ステップと、
前記単結晶(14)内の空孔及び格子間原子の拡散係数及び境界条件に基づいて拡散方程式を解くことにより前記所定の時間間隔の経過した後の空孔及び格子間原子の濃度分布を求める第１０ステップと、
前記空孔の濃度分布に基づいてボイドの形成開始温度を求める第１１ステップと、
前記空孔の濃度分布に基づいて高温酸素析出物の形成開始温度を求める第１２ステップと、
前記単結晶(14)内のそれぞれのメッシュの格子点における温度が次第に低下して前記ボイドの形成開始温度になったときの前記ボイドの濃度を求める第１３ステップと、
前記単結晶(14)内のそれぞれのメッシュの格子点における温度が前記ボイドの形成開始温度より低いときの前記ボイドの半径を求める第１４ステップと、
前記単結晶(14)内のそれぞれのメッシュの格子点における温度が次第に低下して前記高温酸素析出物の形成開始温度になったときの前記高温酸素析出物の濃度を求める第１５ステップと、
前記単結晶(14)内のそれぞれのメッシュの格子点における温度が前記高温酸素析出物の形成開始温度より低いときの前記高温酸素析出物の半径を求める第１６ステップと、
第９ステップから第１６ステップを前記単結晶(14)の冷却が完了するまで繰返す第１７ステップと、
前記高温酸素析出物の存在するＰバンドと前記ボイド及び前記高温酸素析出物のいずれも存在しない無欠陥領域とを区画しかつ前記無欠陥領域が前記Ｐバンドより単結晶トップ側に位置するように区画する第２等濃度線(HC2X)を求める第１８ステップと、
高濃度酸素析出物が存在するＢバンドと前記無欠陥領域とを区画しかつ前記無欠陥領域が前記Ｂバンドより単結晶ボトム側に位置するように区画する第２分布線(BC2X)を求める第１９ステップと、
前記第２等濃度線(HC2X)の変曲点の最小値と前記第２分布線(BC2X)の変曲点の最大値との差(ΔZ2)を求める第２０ステップと、
前記単結晶製造条件のパラメータをＰ2からＰNに順次変え前記第２ステップから前記第２０ステップを繰返して前記第２等濃度線(HC2X)の変曲点の最小値と前記第２分布線(BC2X)の変曲点の最大値との差(ΔZ2)を求めた後に前記差(ΔZ2)が最も大きくなる単結晶製造条件を求める第２１ステップと
を含むコンピュータを用いて単結晶の無欠陥領域を最大化するシミュレーション方法。Parameter groups P 1 , P 2 ,..., P N are arbitrarily defined by varying variables fed back as a single crystal production condition when pulling up the single crystal 14 by the puller 11. The first step;
From the parameter group P 1 , P 2 ,..., P N , the parameter P 1 is used as the single crystal production condition, and the cooling of the single crystal (14) is completed from the start of the pulling of the single crystal (14) from the melt. A second step of modeling the hot zone of the puller (11) up to a mesh structure;
The summarized mesh for each member of the hot zone and the with physical properties of each member for the gathered mesh parameter groups P 1, P 2, ..., first enter the parameters P 1 of the P N to the computer 3 steps,
A fourth step of obtaining the surface temperature distribution of each member based on the heat value of the heater and the radiation rate of each member;
The turbulence obtained by assuming that the melt (12) is turbulent after obtaining the internal temperature distribution of each member by solving the heat conduction equation based on the surface temperature distribution and thermal conductivity of each member. A fifth step of further determining the internal temperature distribution of the melt (12) considering convection by connecting and solving the flow model equation and the Navier-Stokes equations;
A sixth step of obtaining a solid-liquid interface shape of the single crystal (14) and the melt (12) according to an isotherm including a triple point (S) of the single crystal;
The fourth step to the sixth step are repeated until the triple point (S) reaches the melting point of the single crystal (14), and the temperature distribution in the puller (11) is calculated to calculate the single crystal (14). A seventh step of obtaining the coordinates and temperature of the mesh and inputting these data to the computer, respectively;
The pulling length and the pulling height of the single crystal (14) are changed stepwise, and the temperature distribution in the pulling machine (11) is calculated repeatedly from the second step to the seventh step to calculate the single crystal (14 8) to obtain the coordinates and temperature of the mesh and input these data to the computer,
Eighth step of dividing the time from the start of pulling of the single crystal (14) from the melt (12) to the completion of cooling of the single crystal (14) at predetermined intervals. A ninth step of determining the pulling length and height of the single crystal (14) and the temperature distribution in the single crystal (14) from the mesh coordinates and temperature data of the single crystal (14) determined in
The concentration distribution of vacancies and interstitial atoms after the predetermined time interval has been obtained by solving the diffusion equation based on the diffusion coefficients and boundary conditions of vacancies and interstitial atoms in the single crystal (14) 10th step;
An eleventh step of determining a void formation start temperature based on the concentration distribution of the pores;
A twelfth step of determining a formation start temperature of the high-temperature oxygen precipitate based on the concentration distribution of the pores;
A thirteenth step of determining the concentration of the void when the temperature at the lattice point of each mesh in the single crystal (14) gradually decreases to the formation start temperature of the void;
A fourteenth step of determining a radius of the void when a temperature at a lattice point of each mesh in the single crystal (14) is lower than a formation start temperature of the void;
A fifteenth step of determining the concentration of the high-temperature oxygen precipitates when the temperature at the lattice points of each mesh in the single crystal (14) gradually decreases to reach the formation start temperature of the high-temperature oxygen precipitates;
A sixteenth step of obtaining a radius of the high-temperature oxygen precipitate when the temperature at the lattice point of each mesh in the single crystal (14) is lower than the formation start temperature of the high-temperature oxygen precipitate;
A seventeenth step of repeating the ninth to sixteenth steps until the cooling of the single crystal (14) is completed;
The P band where the high-temperature oxygen precipitates are present and the defect-free region where neither the void nor the high-temperature oxygen precipitates exist are partitioned, and the defect-free region is located on the single crystal top side from the P band. An 18th step for obtaining a second isoconcentration line (HC2 X ) to be partitioned;
Obtain a second distribution line (BC2 X ) that partitions the B band in which high-concentration oxygen precipitates are present and the defect-free region and that the defect-free region is positioned on the single crystal bottom side from the B band. 19th step;
A twentieth step of obtaining a difference (ΔZ 2 ) between the minimum value of the inflection point of the second isoconcentration line (HC2 X ) and the maximum value of the inflection point of the second distribution line (BC2 X );
The parameter of the single crystal production condition is sequentially changed from P 2 to PN, and the second step to the twentieth step are repeated, so that the minimum value of the inflection point of the second isoconcentration line (HC2 X ) and the second distribution are obtained. A computer including a twenty-first step for obtaining a single crystal manufacturing condition in which the difference (ΔZ 2 ) is the largest after obtaining the difference (ΔZ 2 ) from the maximum value of the inflection point of the line (BC2 X ). A simulation method that maximizes the defect-free region of a crystal.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2003194779A JP4192704B2 (en) | 2002-07-10 | 2003-07-10 | A simulation method to maximize the defect-free region of a single crystal |
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002200828 | 2002-07-10 | ||
| JP2003194779A JP4192704B2 (en) | 2002-07-10 | 2003-07-10 | A simulation method to maximize the defect-free region of a single crystal |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2004091316A JP2004091316A (en) | 2004-03-25 |
| JP4192704B2 true JP4192704B2 (en) | 2008-12-10 |
Family
ID=32072138
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2003194779A Expired - Lifetime JP4192704B2 (en) | 2002-07-10 | 2003-07-10 | A simulation method to maximize the defect-free region of a single crystal |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP4192704B2 (en) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP6135611B2 (en) * | 2014-07-03 | 2017-05-31 | 信越半導体株式会社 | Point defect concentration calculation method, Grown-in defect calculation method, Grown-in defect in-plane distribution calculation method, and silicon single crystal manufacturing method using them |
| CN113744818B (en) * | 2021-09-03 | 2023-09-15 | 上海大学 | Prediction method for ternary rare earth oxide composite point defects |
-
2003
- 2003-07-10 JP JP2003194779A patent/JP4192704B2/en not_active Expired - Lifetime
Also Published As
| Publication number | Publication date |
|---|---|
| JP2004091316A (en) | 2004-03-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US6045610A (en) | Methods of manufacturing monocrystalline silicon ingots and wafers by controlling pull rate profiles in a hot zone furnance | |
| EP1310583B1 (en) | Method for manufacturing of silicon single crystal wafer | |
| JP4380537B2 (en) | Method for producing silicon single crystal | |
| WO1994016124A1 (en) | Method and apparatus for predicting crystal quality of single-crystal semiconductor | |
| JPH11199387A (en) | Production of silicon single crystal and silicon single crystal wafer | |
| JP4701738B2 (en) | Single crystal pulling method | |
| JP4548306B2 (en) | Method for producing silicon single crystal | |
| US20120279438A1 (en) | Methods for producing single crystal silicon ingots with reduced incidence of dislocations | |
| JP4192704B2 (en) | A simulation method to maximize the defect-free region of a single crystal | |
| JP6052189B2 (en) | Heat treatment method for silicon single crystal wafer | |
| JP2003055092A (en) | Method of pulling silicon single crystal | |
| EP1591566B1 (en) | Method of producing p-doped silicon single crystal and p-doped n-type silicon single crystal wafe | |
| JP4154936B2 (en) | Single crystal defect-free region simulation method | |
| KR100665683B1 (en) | Method for manufacturing silicon single crystal | |
| JP4102966B2 (en) | Pulling method of silicon single crystal | |
| JP3752890B2 (en) | Method for producing silicon single crystal ingot | |
| JP2001089294A (en) | Method of continuously pulling up silicon single crystal free from agglomerates of point defects | |
| Virzi et al. | Computer simulation of Czochralski silicon thermal history and its effect on bulk stacking fault nuclei generation | |
| JP4449347B2 (en) | OSF ring distribution prediction method by simulation | |
| JP4182691B2 (en) | A method for predicting pulling conditions of pure silicon single crystals. | |
| WO2004106594A1 (en) | Method of simulation with respect to density distribution and size distribution of void defect within single crystal and oxygen precipitation nucleus within single crystal | |
| JP2003055091A (en) | Method of pulling silicon single crystal | |
| JP4403722B2 (en) | Simulation method for density distribution and size distribution of void defect in silicon single crystal | |
| JP2004363412A (en) | Method for simulating density distribution and size distribution of oxygen precipitation nuclei in single crystal | |
| JP4193470B2 (en) | A method for predicting the relationship between the flow rate of inert gas and the pure margin. |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A625 | Written request for application examination (by other person) |
Free format text: JAPANESE INTERMEDIATE CODE: A625 Effective date: 20050704 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20070730 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20080826 |
|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20080908 |
|
| R150 | Certificate of patent or registration of utility model |
Ref document number: 4192704 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111003 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20121003 Year of fee payment: 4 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20131003 Year of fee payment: 5 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| EXPY | Cancellation because of completion of term |