JP4604462B2

JP4604462B2 - Simulation method of density distribution and size distribution of oxygen precipitation nuclei in single crystal

Info

Publication number: JP4604462B2
Application number: JP2003161493A
Authority: JP
Inventors: 浩之介北村; 純古川; 直樹小野
Original assignee: Sumco Corp
Current assignee: Sumco Corp
Priority date: 2003-06-06
Filing date: 2003-06-06
Publication date: 2011-01-05
Anticipated expiration: 2023-06-06
Also published as: JP2004363412A

Description

【０００１】
【発明の属する技術分野】
本発明は、チョクラルスキー（以下、ＣＺという。）法にて引上げられるシリコン等の単結晶内の酸素析出核の密度分布及びサイズ分布をコンピュータシミュレーションする方法に関するものである。ここで言う酸素析出核とは、単結晶内の空孔と酸素を消費して形成された酸素析出核であり、これを低温（７００〜８００℃）の熱処理で安定化し、更に高温（１０００〜１１００℃）の熱処理で成長させると、光学顕微鏡などでＢＭＤ（Bulk Micro-Defect）として観察される。
【０００２】
【従来の技術】
従来、この種のシミュレーション方法として、図７に示すように、総合伝熱シミュレータを用いてＣＺ法によるシリコン単結晶４引上げ時の引上げ機１内のホットゾーン構造及びそのシリコン単結晶４の引上げ速度に基づいて、シリコン融液２の熱伝導率を操作することによりシリコン融液２の内部温度分布を予測し、この内部温度分布からシリコン単結晶４のメッシュの座標及び温度をそれぞれ求め、更にシリコン単結晶４内の格子間シリコン及び空孔の拡散係数及び境界条件に基づいて拡散方程式を解くことにより、上記格子間シリコン及び空孔の密度分布をコンピュータを用いて求める方法が知られている。このシミュレーション方法では、ホットゾーンの各部材がメッシュ分割されてモデル化される。特にシリコン融液２のメッシュは計算時間を短くするために１０ｍｍ程度と比較的粗く設定される。
【０００３】
一方、ＣＺ法によるシリコン単結晶育成時の操業条件及び炉内の温度条件をコンピュータにデータとして入力し、単結晶内部に導入される空孔起因のグローンイン欠陥のサイズ及び密度と単結晶の成長速度との関係をコンピュータを用いたシミュレーションにより求め、このシミュレーション結果に基づいて単結晶内部に導入される空孔起因のグローンイン欠陥のサイズ及び密度が所定値になるように成長速度を選択して単結晶を育成するシリコン単結晶の製造方法が知られている（例えば、特許文献１参照）。このシリコン単結晶の製造方法では、総合伝熱解析プログラムを用いて炉内の熱分布のシミュレーションを行い、単結晶の育成時の単結晶中心の鉛直方向における温度勾配を求める。
【０００４】
【特許文献１】
特開２００２−１４５６９６号公報
【０００５】
【発明が解決しようとする課題】
しかし、上記従来の格子間シリコン及び空孔の密度分布のシミュレーション方法では、実際の引上げ機においては発生するシリコン融液の対流を考慮しておらず、またシリコン融液のメッシュが比較的粗いため、固液界面形状の再現性が悪く、精度の良い単結晶内温度分布を提供できない。そのため、シリコン単結晶内の格子間シリコン及び空孔の密度分布が実測値と大幅に相違してしまい、結晶冷却過程に比較的低温で生成されるシリコン単結晶内の欠陥である酸素析出核の密度分布及びサイズ分布は判らなかった。
一方、上記特許文献１に示されたシリコン単結晶の製造方法では、シリコン単結晶の鉛直方向のグローンイン欠陥のサイズ分布及び密度分布を予測できるけれども、シリコン単結晶の半径方向のグローイン欠陥のサイズ分布及び密度分布を予測できず、ウェーハ内面のグローイン欠陥を検討できない問題点があった。
【０００６】
本発明の目的は、融液の対流を考慮して成長中の単結晶内の温度分布を解析した後に、融液から切離された単結晶の冷却過程を考慮して解析することにより、単結晶内の酸素析出核の密度分布及びサイズ分布を正確に予測できる、単結晶内酸素析出核の密度分布及びサイズ分布のシミュレーション方法を提供することにある。
【０００７】
【課題を解決するための手段】
請求項１に係る発明は、図１〜図６に示すように、引上げ機１１による単結晶１４の融液１２からの引上げ開始時から単結晶１４の冷却完了時までの引上げ機１１のホットゾーンをメッシュ構造でモデル化する第１ステップと、ホットゾーンの各部材毎にメッシュをまとめかつこのまとめられたメッシュに対する各部材の物性値とともに単結晶１４の引上げ長及びこの引上げ長に対応する単結晶１４の引上げ速度をそれぞれコンピュータに入力する第２ステップと、各部材の表面温度分布をヒータの発熱量及び各部材の輻射率に基づいて求める第３ステップと、各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式を解くことにより各部材の内部温度分布を求めた後に融液１２が乱流であると仮定して得られた乱流モデル式及びナビエ・ストークスの方程式を連結して解くことにより対流を考慮した融液１２の内部温度分布を更に求める第４ステップと、単結晶１４及び融液１２の固液界面形状を単結晶の三重点Ｓを含む等温線に合せて求める第５ステップと、第３ステップから第５ステップを三重点Ｓが単結晶１４の融点になるまで繰返し引上げ機１１内の温度分布を計算して単結晶１４のメッシュの座標及び温度を求めこれらのデータをそれぞれコンピュータに入力する第６ステップと、単結晶１４の引上げ長及び引上げ高さを段階的に変えて第１ステップから第６ステップまでを繰返し引上げ機１１内の温度分布を計算して単結晶１４のメッシュの座標及び温度を求めこれらのデータをそれぞれコンピュータに入力する第７ステップと、単結晶１４の融液１２からの引上げ開始時から単結晶１４の冷却完了時までの時間を所定の間隔毎に区切りこの区切られた時間間隔毎に第７ステップで求めた単結晶１４のメッシュの座標及び温度のデータから単結晶１４の引上げ長及び引上げ高さと単結晶１４内の温度分布とを求める第８ステップと、単結晶１４内の空孔及び格子間原子の拡散係数及び境界条件に基づいて拡散方程式を解くことにより所定の時間間隔の経過した後の空孔及び格子間原子の密度分布を求める第９ステップと、空孔の密度分布に基づいてボイド２１の形成開始温度を求める第１０ステップと、単結晶１４内のそれぞれのメッシュの格子点における温度が次第に低下してボイド２１の形成開始温度になったときのボイド２１の密度を求める第１１ステップと、単結晶１４内のそれぞれのメッシュの格子点における温度がボイド２１の形成開始温度より低いときのボイド２１の半径を求める第１２ステップと、ボイド２１の周囲に成長する内壁酸化膜２２の形成開始温度を求める第１３ステップと、単結晶１４内のそれぞれのメッシュの格子点における温度が次第に低下して内壁酸化膜２２の形成開始温度より低くなったときのボイド２１の半径と内壁酸化膜２２の厚さとを互いに関連させて求める第１４ステップと、単結晶１４内の酸素の密度分布と空孔の密度分布と温度分布に基づいて酸素析出核の単位時間当りの発生量と臨界半径とを求める第１５ステップと、酸素析出核の単位時間当りの発生量がゼロを越えたときにボイド２１の半径及び内壁酸化膜２２の厚さを求める計算を停止しかつ酸素析出核の半径を求める第１６ステップと、第８ステップから第１６ステップを単結晶１４の冷却が完了するまで繰返す第１７ステップとを含む、コンピュータを用いて単結晶内酸素析出核の密度分布及びサイズ分布のシミュレーションを行う方法である。
【０００８】
この請求項１に記載されたシミュレーション方法では、融液１２の対流を考慮して融液１２から成長する単結晶１４内の温度分布を求めるだけでなく、更に冷却過程における単結晶１４内の温度分布までも求めることによって、単結晶１４内の酸素析出核の密度分布及びサイズ分布を正確に予測できる。即ち、融液１２から切離された単結晶１４の冷却過程における単結晶１４の徐冷及び急冷の効果を考慮して解析することによって、ボイド２１及び内壁酸化膜２２からなるボイド欠陥の密度分布とサイズ分布を求め、酸素析出核の単位時間当りの発生量とこの臨界半径を求め、更に酸素析出核の単位時間当りの発生量がゼロを越えたときに、ボイド半径及び内壁酸化膜厚さを求める計算を停止し、かつ酸素析出核の半径を求める。これにより単結晶１４内の酸素析出核の密度分布及びサイズ分布を正確に予測できる。
ここで、酸素析出核の臨界半径とは、結晶成長中に発生した核の最小半径をいい、この半径は空孔密度の過飽和度及び酸素密度の過飽和度に依存する。
【０００９】
【発明の実施の形態】
次に本発明の実施の形態を図面に基づいて説明する。
図５に示すように、シリコン単結晶引上げ機１１のチャンバ内には、シリコン融液１２を貯留する石英るつぼ１３が設けられる。この石英るつぼ１３は図示しないが黒鉛サセプタ及び支軸を介してるつぼ駆動手段に接続され、るつぼ駆動手段は石英るつぼ１３を回転させるとともに昇降させるように構成される。また石英るつぼ１３の外周面は石英るつぼ１３から所定の間隔をあけてヒータ（図示せず）により包囲され、このヒータは保温筒（図示せず）により包囲される。ヒータは石英るつぼ１３に投入された高純度のシリコン多結晶体を加熱・溶融してシリコン融液１２にする。またチャンバの上端には図示しないが円筒状のケーシングが接続され、このケーシングには引上げ手段が設けられる。引上げ手段はシリコン単結晶１４を回転させながら引上げるように構成される。
【００１０】
このように構成されたシリコン単結晶引上げ機１１におけるシリコン単結晶１４内の酸素析出核の密度分布及びサイズ分布のシミュレーション方法を図１〜図６に基づいて説明する。
先ず第１ステップとして、シリコン単結晶１４を所定長さＬ₁（例えば１００ｍｍ）まで引上げた状態におけるシリコン単結晶引上げ機１１のホットゾーンの各部材、即ちチャンバ，石英るつぼ１３，シリコン融液１２，シリコン単結晶１４，黒鉛サセプタ，保温筒等をメッシュ分割してモデル化する。具体的には上記ホットゾーンの各部材のメッシュ点の座標データをコンピュータに入力する。このときシリコン融液１２のメッシュのうちシリコン単結晶１４の径方向のメッシュであってかつシリコン融液１２のシリコン単結晶１４直下の一部又は全部のメッシュ（以下、径方向メッシュという。）を０．０１〜５．００ｍｍ、好ましくは０．２５〜１．００ｍｍに設定する。またシリコン融液１２のメッシュのうちシリコン単結晶１４の長手方向のメッシュであってかつシリコン融液１２の一部又は全部のメッシュ（以下、長手方向メッシュという。）を０．０１〜５．００ｍｍ、好ましくは０．１〜０．５ｍｍに設定する。
【００１１】
径方向メッシュを０．０１〜５．００ｍｍの範囲に限定したのは、０．０１ｍｍ未満では計算時間が極めて長くなり、５．００ｍｍを越えると計算が不安定になり、繰返し計算を行っても固液界面形状が一定に定まらなくなるからである。
また長手方向メッシュを０．０１〜５．００ｍｍの範囲に限定したのは、０．０１ｍｍ未満では計算時間が極めて長くなり、５．００ｍｍを越えると固液界面形状の計算値が実測値と一致しなくなるからである。なお、径方向メッシュの一部を０．０１〜５．００の範囲に限定する場合には、シリコン単結晶１４直下のシリコン融液１２のうちシリコン単結晶１４外周縁近傍のシリコン融液１２を上記範囲に限定することが好ましく、長手方向メッシュの一部を０．０１〜５．００の範囲に限定する場合には、シリコン融液１２の液面近傍及び底近傍を上記範囲に限定することが好ましい。
【００１２】
第２ステップとして上記ホットゾーンの各部材毎にメッシュをまとめ、かつこのまとめられたメッシュに対して各部材の物性値をそれぞれコンピュータに入力する。例えば、チャンバがステンレス鋼にて形成されていれば、そのステンレス鋼の熱伝導率，輻射率，粘性率，体積膨張係数，密度及び比熱がコンピュータに入力される。またシリコン単結晶１４の引上げ長及びこの引上げ長に対応するシリコン単結晶１４の引上げ速度と、後述する乱流モデル式（１）の乱流パラメータＣとをコンピュータに入力する。
【００１３】
第３ステップとして、ホットゾーンの各部材の表面温度分布をヒータの発熱量及び各部材の輻射率に基づいてコンピュータを用いて求める。即ち、ヒータの発熱量を任意に設定してコンピュータに入力するとともに、各部材の輻射率から各部材の表面温度分布をコンピュータを用いて求める。次に第４ステップとしてホットゾーンの各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式（１）をコンピュータを用いて解くことにより各部材の内部温度分布を求める。ここでは、記述を簡単にするためｘｙｚ直交座標系を用いたが、実際の計算では円筒座標系を用いる。
【数１】

【００１４】
上記式（１）において、ρは各部材の密度、ｃは各部材の比熱、Ｔは各部材の各メッシュ点での絶対温度、ｔは時間である。またλ_x，λ_y及びλ_zは各部材の熱伝導率のｘ，ｙ及びｚ方向成分であり、ｑはヒータの発熱量である。
【００１５】
一方、シリコン融液１２に関しては、上記熱伝導方程式（１）でシリコン融液１２の内部温度分布を求めた後に、このシリコン融液１２の内部温度分布に基づき、シリコン融液１２が乱流であると仮定して得られた乱流モデル式（２）及びナビエ・ストークスの方程式（３）〜（５）を連結して、シリコン融液１２の内部流速分布をコンピュータを用いて求める。
【数２】

【００１６】
上記式（２）において、κ_tはシリコン融液１２の乱流熱伝導率であり、ｃはシリコン融液１２の比熱であり、Ｐｒ_tはプラントル数であり、ρはシリコン融液１２の密度であり、Ｃは乱流パラメータであり、ｄはシリコン融液１２を貯留する石英るつぼ１３壁からの距離であり、ｋはシリコン融液１２の平均流速に対する変動成分の二乗和である。
【数３】

【００１７】
上記式（３）〜式（５）において、ｕ，ｖ及びｗはシリコン融液１２の各メッシュ点での流速のｘ，ｙ及びｚ方向成分であり、ν_lはシリコン融液１２の分子動粘性係数（物性値）であり、ν_tはシリコン融液１２の乱流の効果による動粘性係数であり、Ｆ_x，Ｆ_y及びＦ_zはシリコン融液１２に作用する体積力のｘ，ｙ及びｚ方向成分である。
【００１８】
上記乱流モデル式（２）はｋｌ（ケイエル）−モデル式と呼ばれ、このモデル式の乱流パラメータＣは０．４〜０．６の範囲内の任意の値が用いられることが好ましい。乱流パラメータＣを０．４〜０．６の範囲に限定したのは、０．４未満又は０．６を越えると計算により求めた界面形状が実測値と一致しないという不具合があるからである。また上記ナビエ・ストークスの方程式（３）〜（５）はシリコン融液１２が非圧縮性であって粘度が一定である流体としたときの運動方程式である。
【００１９】
上記求められたシリコン融液１２の内部流速分布に基づいて熱エネルギー方程式（６）を解くことにより、シリコン融液１２の対流を考慮したシリコン融液１２の内部温度分布をコンピュータを用いて更に求める。
【数４】

【００２０】
上記式（６）において、ｕ，ｖ及びｗはシリコン融液１２の各メッシュ点での流速のｘ，ｙ及びｚ方向成分であり、Ｔはシリコン融液１２の各メッシュ点での絶対温度であり、ρはシリコン融液１２の密度であり、ｃはシリコン融液１２の比熱であり、κ_lは分子熱伝導率（物性値）であり、κ_tは式（１）を用いて計算される乱流熱伝導率である。
【００２１】
次いで第５ステップとして、シリコン単結晶１４及びシリコン融液１２の固液界面形状を図５の点Ｓで示すシリコンの三重点Ｓ（固体と液体と気体の三重点(tri-junction)）を含む等温線に合せてコンピュータを用いて求める。第６ステップとして、コンピュータに入力するヒータの発熱量を変更し（次第に増大し）、上記第３ステップから第５ステップを三重点がシリコン単結晶１４の融点になるまで繰返した後に、引上げ機１１内の温度分布を計算してシリコン単結晶のメッシュの座標及び温度を求め、これらのデータをコンピュータに記憶させる。
【００２２】
次に第７ステップとして、シリコン単結晶１４の引上げ長Ｌ₁にδ（例えば５０ｍｍ）だけ加えて上記第１ステップから第６ステップまでを繰返した後に、引上げ機１１内の温度分布を計算してシリコン単結晶１４のメッシュの座標及び温度を求め、これらのデータをコンピュータに記憶させる。この第７ステップは、シリコン単結晶１４の引上げ長Ｌ₁が長さＬ₂（Ｌ₂はシリコン融液１２から切離されたときのシリコン単結晶１４の長さ（成長完了時の結晶長）である。）に達してシリコン単結晶１４がシリコン融液１２から切離された後、更にシリコン単結晶１４が引上げられてその高さＨ₁（Ｈ₁はシリコン単結晶１４の直胴開始部からシリコン融液１２の液面までの距離である（図５）。）がＨ₂（Ｈ₂は冷却完了時のシリコン単結晶１４の直胴開始部からシリコン融液１２の液面までの距離である。）に達するまで、即ちシリコン単結晶１４の冷却が完了するまで行われる。なお、シリコン単結晶１４がシリコン融液１２から切離された後は、シリコン単結晶１４の引上げ高さＨ₁にδ（例えば５０ｍｍ）だけ加え、上記と同様に上記第１ステップから第６ステップまでを繰返す。
【００２３】
シリコン単結晶１４の引上げ高さＨ₁がＨ₂に達すると、第８ステップに移行する。第８ステップでは、シリコン単結晶１４をシリコン融液１２から成長させて引上げ始めたときｔ₀から、シリコン単結晶１４をシリコン融液１２から切離して更にシリコン単結晶１４を引上げ、その冷却が完了したときｔ₁までの時間を、所定の間隔Δｔ秒（微小時間間隔）毎に区切る。このときシリコン単結晶１４内の格子間シリコン及び空孔の拡散係数及び境界条件のみならず、後述するボイド２１及び内壁酸化膜２２（図６）の密度分布及びサイズ分布を求めるための式に用いられる定数をそれぞれコンピュータに入力する。上記区切られた時間間隔Δｔ秒毎に、第７ステップで求めたシリコン単結晶１４のメッシュの座標及び温度のデータから、シリコン単結晶１４の引上げ長Ｌ₁及び引上げ高さＨ₁と、シリコン単結晶１４内の温度分布とを求める。
【００２４】
即ち、第１〜第７ステップでシリコン単結晶のメッシュの座標及び温度を引上げ長δ毎に求め、シリコン単結晶を例えば５０ｍｍ引上げるのに数十分要するため、この数十分間でのシリコン単結晶のメッシュの温度変化を時間の関数として微分することにより、時刻ｔ₀からΔｔ秒後におけるシリコン単結晶１４の引上げ長Ｌ₁及び引上げ高さＨ₁とシリコン単結晶１４内の温度分布を算出する。次にシリコン単結晶１４内の空孔及び格子間シリコンの拡散係数及び境界条件に基づいて拡散方程式を解くことにより、Δｔ秒経過後の空孔及び格子間シリコンの密度分布を求める（第９ステップ）。
【００２５】
具体的には、空孔の密度Ｃ_vの計算式が次の式（７）で示され、格子間シリコンの密度Ｃ_iの計算式が次の式（８）で示される。式（７）及び式（８）において、密度Ｃ_v及び密度Ｃ_iの経時的進展を計算するために、空孔と格子間シリコンの熱平衡がシリコン単結晶の全表面で維持されると仮定する。
【数５】

【００２６】
上記式（７）及び式（８）において、Ｋ₁及びＫ₂は定数であり、Ｅ_i及びＥ_vはそれぞれ格子間シリコン及び空孔の形成エネルギーであり、Ｃ_ve及びＣ_ieは空孔の平衡密度及び格子間シリコンの平衡密度である。またｋ_Bはボルツマン定数、Ｔは絶対温度を意味する。
【００２７】
上記平衡式は時間で微分され、空孔及び格子間シリコンに対してそれぞれ次の式（９）及び式（１０）になる。
【数６】

【００２８】
上記式（９）及び式（１０）において、Θ(ｘ)はヘビサイド関数（Heaviside function）である。即ち、ｘ＜０のときΘ(ｘ)＝０であり、かつｘ＞０のときΘ(ｘ)＝１である。またＴ_pは内壁酸化膜２２（ボイド２１表面に形成されたシリコンの酸化膜（ＳｉＯｘ膜））の形成開始温度であり、Ｔ_vはボイド２１の形成開始温度である（図６）。Ｄ_Oは酸素拡散係数であり、Ｃ_Oは酸素析出核形成温度での酸素密度である。更にＪ_cは臨界半径を有する酸素析出核の単位時間当りの発生量であり、ｒ_cは酸素析出核の半径である。なお、式（９）及び（１０）のそれぞれ右側の第１項はフィックの拡散式であり、右側の第１項中のＤ_v及びＤ_iは、次の式（１１）及び（１２）で表される空孔及び格子間シリコンの拡散係数である。
【数７】

【００２９】
上記式（１１）及び式（１２）において、△Ｅ_v及び△Ｅ_iはそれぞれ空孔及び格子間シリコンの活性化エネルギーであり、ｄ_v及びｄ_iはそれぞれ定数である。また式（９）及び式（１０）のそれぞれ右側の第２項中のＥ_v ^t及びＥ_i ^tは熱拡散による空孔及び格子間シリコンの活性化エネルギーであり、ｋ_Bはボルツマン定数である。式（９）及び式（１０）のそれぞれ右側の第３項のｋ_ivは空孔及び格子間シリコンペアの再結合定数である。式（９）及び式（１０）のそれぞれ右側の第４項のＮ_vはボイド２１の密度であり、ｒ_vはボイド２１の半径であり、更に式（９）の右側の第５項のＮ_pは内壁酸化膜２２の密度であり、Ｒ_pは内壁酸化膜２２の外半径である。
【００３０】
上記式（９）が成立つのは、空孔が析出するための空孔の流束が十分に大きく、シリコン単結晶１４を構成するＳｉマトリックスと内壁酸化膜を構成するＳｉＯｘとの単位質量当りの体積差を埋められる場合、即ちγＤ_OＣ_O≦Ｄ_v（Ｃ_v−Ｃ_ve）の場合である。上記以外の場合、即ちγＤ_OＣ_O＞Ｄ_v（Ｃ_v−Ｃ_ve）の場合には、次式（１３）が成立つ。ここで、γは酸化物が歪みなく成長するのに必要な酸素１原子当たりの空孔量であり、酸素析出核の相に依存する数である。酸素析出核がＳｉＯ₂からなる場合にはγ≒０．５であり、酸素析出核がＳｉＯからなる場合にはγ≒０．６５である。またＤ_Oは酸素の拡散係数であり、Ｃ_Oは酸素析出核形成温度での酸素密度である。
【数８】

【００３１】
上記式（１３）において、Θ(ｘ)はヘビサイド関数（Heaviside function）である。即ち、ｘ＜０のときΘ(ｘ)＝０であり、かつｘ＞０のときΘ(ｘ)＝１である。またＴ_pは内壁酸化膜２２の形成開始温度であり、Ｔ_vはボイド２１の形成開始温度である。更にＪ_cは臨界半径を有する酸素析出核の単位時間当りの発生量であり、ｒ_cは酸素析出核の半径である。なお、式（１３）の右側の第１項はフィックの拡散式であり、右側の第１項中のＤ_vは、式（１１）で表される空孔の拡散係数である。
【００３２】
次に第１０ステップとして、上記拡散方程式を解くことにより求めた空孔の密度Ｃ_v分布に基づいて、ボイド２１の形成開始温度Ｔ_vを次の式（１４）から求める。
【数９】

【００３３】
上記式（１４）において、Ｃ_vはシリコン単結晶１４中の空孔密度であり、Ｃ_v0はシリコン単結晶１４の融点Ｔ_mでの空孔平衡密度であり、Ｅ_vは空孔形成エネルギーである。またσ_vはシリコン単結晶１４の結晶面(１１１)における界面エネルギーであり、ρはシリコン単結晶１４の密度であり、ｋ_Bはボルツマン定数である。
【００３４】
第１１ステップとして、シリコン単結晶１４内のそれぞれのメッシュの格子点における温度が次第に低下してボイド２１の形成開始温度Ｔ_vになったときに、次の近似式（１５）を用いてボイド２１の密度Ｎ_vを求める。
【数１０】

【００３５】
上記式（１５）において、ｎ_vはボイド２１の密度であり、（dT/dt）はシリコン単結晶１４の冷却速度である。上記ボイド２１の密度は空孔過飽和度に依存するけれども、計算するシリコン単結晶１４での核形成温度のような狭い温度領域では一定とみなしてよい。またＤ_vはボイド２１の拡散係数であり、ｋ_Bはボルツマン定数であり、Ｃ_vはシリコン単結晶１４中の空孔密度である。
【００３６】
第１２ステップとして、シリコン単結晶１４内のそれぞれのメッシュの格子点における温度がボイド２１の形成開始温度Ｔ_vより低いときのボイドの半径ｒ_vを、次の式（１６）から求める。ここでボイド２１の成長速度は空孔の拡散速度に依存する拡散律速である。またボイド２１の形状は実際には八面体であるけれども、ここでは計算の効率化から球状として扱う。
【数１１】

【００３７】
上記式（１６）において、ｔ１はシリコン単結晶１４のメッシュの格子点における温度がボイド２１の形成開始温度Ｔ_vまで低下したときの時刻である。またｒ_vcrはボイド２１の臨界半径であり、この臨界半径ｒ_vcrは上記時刻ｔ１での値とする。更にＣ_ve及びＣ_ieはそれぞれの状態での空孔及び格子間シリコンの平衡密度である。
【００３８】
第１３ステップとして、ボイド２１の周囲に成長する内壁酸化膜２２の形成開始温度Ｔ_pを求める。内壁酸化膜２２の形成開始温度Ｔ_pがボイド２１の形成開始温度Ｔ_vより小さい場合には、酸素と空孔が結合し、ボイド２１表面に内壁酸化膜２２が成長する（図６）。本発明のモデルでは、上記内壁酸化膜２２はボイド２１が発生するとすぐに成長するものとして扱った。ここで、「すぐに」とは、「ボイドが発生した次の繰返し計算から」という意味である。
【００３９】
第１４ステップとして、シリコン単結晶１４内のそれぞれのメッシュの格子点における温度が次第に低下して内壁酸化膜２２の形成開始温度Ｔ_pより低くなったときのボイド２１の半径ｒ_vと内壁酸化膜２２の厚さｄとを互いに関連させて求める。一般的に、酸素はその拡散とシリコンとの反応によって内壁酸化膜２２が形成される。酸素が内壁酸化膜２２外からこの酸化膜の外周面に流入する流束Ｊ_Oは、ボイド２１の中心を原点とする内壁酸化膜２２の外半径をＲ_pとするとき、次の式（１７）から求められる。
【数１２】

【００４０】
上記式（１７）において、Ｄ_Oは酸素の拡散係数であり、Ｃ_Oは酸素密度であり、Ｃ_Oeifは球状のボイド２１の中心を原点とする半径Ｒ_pにおける酸素平衡密度であり、ｇ_OはＳｉＯｘ及びＳｉの界面、即ち内壁酸化膜２２の外周面における酸素原子の反応割合である。また空孔が内壁酸化膜２２外からこの酸化膜の外周面（半径Ｒ_pの球面）に流入する流束Ｊ_vは次の式（１８）から求められる。
【数１３】

【００４１】
上記式（１８）において、ｇ_vは内壁酸化膜２２の外周面での空孔の反応割合であり、Ｃ_veifは内壁酸化膜２２の外周面での空孔の平衡密度である。シリコンの内壁酸化膜２２が成長すると同時に空孔が消費され、空孔の消費割合は酸素１原子に対してγである。このため内壁酸化膜２２が歪みなく成長するための条件は、酸素の流束Ｊ_Oに対して空孔の流束はγＪ_Oであるので、空孔の流束γＪ_Oは次の式（１９）から求められる。
【数１４】

【００４２】
従って、(Ｊ_v−γＪ_O)は内壁酸化膜２２に吸収される空孔の流束となる。ここで空孔の内壁酸化膜２２外からこの酸化膜への流束が増大する場合((Ｊ_v−γＪ_O)＞０)、或いは空孔のボイド２１内から酸化膜への流束が時間の経過に対して増大する場合((Ｊ_v−γＪ_O)＜０)の上記流束の時間に対する変化の割合をαとすると、このαは内壁酸化膜２２の厚さに依存すると考えられる。このため内壁酸化膜２２の厚さｄ、即ち(Ｒ_p−ｒ_v)がｄ₀未満であるならば、内壁酸化膜２２は空孔を酸化膜外及びボイド２１内の双方向から部分浸透させると考えられる。換言すれば、ｄ＜ｄ₀であってα≠０である場合には、空孔はある割合で内壁酸化膜２２を貫通してボイド２１の成長に使われる。但し、格子間シリコンは酸素析出層が核形成されると同時にボイド２１の成長に使われなくなる。
【００４３】
また空孔の内壁酸化膜２２外からこの酸化膜への流束量が少ない場合、即ち(Ｊ_v−γＪ_O)＜０である場合には、歪みの無い内壁酸化膜２２を成長させることができないので、空孔はボイド２１内から酸化膜に流入し、ボイド２１の外径は小さくなる。この結果、ｄ≧ｄ₀であってα＝０の場合には、内壁酸化膜２２は空孔を全く通過させなくなる。そして、n_vをボイド２１中の空孔の数とすると、dn_v/dtは空孔の流入量となり、dn_v/dt＝α(Ｊ_v−γＪ_O)であるならば、次の式（２０）が成り立つ。なお、式（２０）において、ｒ_vは式（１６）で求めたボイド２１の半径であり、Ｒ_pは内壁酸化膜２２の外半径である。またＤ_Oは酸素の拡散係数であり、Ｃ_Oは酸素密度である。更にＣ_Oeifは球状のボイド２１の中心を原点とする半径Ｒ_pにおける酸素平衡密度であり、ｇ_Oは内壁酸化膜２２の外周面における酸素原子の反応割合である。
【数１５】

【００４４】
同時に内壁酸化膜２２を構成するＳｉＯｘ分子の数は式（２０）の酸素が内壁酸化膜２２外からこの酸化膜の外周面に流入する流束Ｊ_Oの消費に依存して変化するため、次の式（２１）が成り立つ。なお、式（２１）において、ρ_p＝ｘ／Ω_SiOxと定義する。このρ_pの定義式において、Ω_SiOxはＳｉＯｘの分子容量（molecular volume）であり、その単位は分子量／密度である。またｘはＳｉＯｘのｘに対応する。
【数１６】

【００４５】
上記式（２１）を上記式（２０）に代入すると次の式（２２）が得られる。なお、式（２２）において、ηはΩ_SiOx／Ω_Siである。このηの定義式において、Ω_Siは１／ρであり、ρはシリコン単結晶１４の密度である。
【数１７】

【００４６】
上記式（２０）及び式（２２）から次のことが分かる。先ずｄがｄ₀になった時刻をｔ₃とすると、この時刻ｔ₃においてαがゼロとなるため、内壁酸化膜２２は非浸透となる。またボイド２１の外径がｒ_v(ｔ₃)で固定されている間は、内壁酸化膜２２の成長はシリコン単結晶１４のマトリックスから流入する空孔或いは酸素の消費量によって決まる。即ち、（Ｊ_v−γＪ_O）＞０である場合には式（２２）は次の式（２３）となり、（Ｊ_v−γＪ_O）＜０である場合には式（２２）は次の式（２４）となる。
【数１８】

【００４７】
上記式（２３）及び式（２４）において、Ω_SiOxはＳｉＯｘの分子容量（molecular volume）である。また式（２３）及び式（２４）において、内壁酸化膜の外周面での空孔平衡密度はＲ_pと無関係であり、Ｃ_veif＝Ｃ_veかつＣ_Oeif＝Ｃ_Oeであると想定している。式（２３）又は式（２４）から内壁酸化膜２２の外半径Ｒ_pを求め、内壁酸化膜２２の外半径Ｒ_pからボイド２１の半径ｒ_vを引くことにより、内壁酸化膜２２の厚さｄを求める。
【００４８】
次に第１５ステップとして、シリコン単結晶１４内の酸素の密度分布と空孔の密度分布と温度分布に基づいて酸素析出核の単位時間当りの発生量Ｊ_cと上記臨界半径ｒ_crとを求める。具体的には、酸素析出核成長（ＳｉＯｘ相）は、酸素密度Ｃ_Oと空孔密度Ｃ_vと温度Ｔに依存するため、古典論では定常状態での酸素析出核の単位時間当りの発生量Ｊ_cを次の式（２５）から求める。
【数１９】

【００４９】
上記式（２５）において、σ_cは酸素析出核形成温度でのＳｉＯｘとＳｉとの界面エネルギーであり、Ｄ_Oは酸素拡散係数である。またＣ_Oは酸素析出核形成温度での酸素密度であり、ｋ_Bはボルツマン定数である。更にＦ^*は酸素析出核の成長バリアであり、酸素析出核を球体と仮定すれば、Ｆ^*は次の式（２６）により求めることができる。
【数２０】

【００５０】
上記式（２６）において、ρ_cはＳｉＯｘからなる酸素析出核の酸素数密度であり、ｆ_cは酸素析出核形成のために必要な駆動力（酸素原子１個当りの力）であり、このｆ_cを次の式（２７）から求める。
【数２１】

【００５１】
上記式（２７）において、Ｃ_Oは酸素析出核形成温度での酸素密度であり、Ｃ_Oeは酸素平衡密度である。またγは酸化物が歪みなく成長するのに必要な酸素１原子当たりの空孔量であり、酸素析出核の相に依存する数である。酸素析出核がＳｉＯ₂からなる場合にはγ≒０．５であり、酸素析出核がＳｉＯからなる場合にはγ≒０．６５である。更にＣ_vは空孔の密度であり、Ｃ_veは空孔の平衡密度である。なお、式（２５）及び式（２６）におけるσ_c（酸素析出核形成温度でのＳｉＯｘとＳｉとの界面エネルギー）は核形成を支配する重要なパラメータであり、このパラメータσ_cは核形成に依存しており、一般的には温度と酸素析出核の半径によって値が変化する。一方、酸素析出核の臨界半径ｒ_crを次の式（２８）により求める。
【数２２】

【００５２】
次に第１６ステップとして、酸素析出核の単位時間当りの発生量がゼロを越えると、ボイド２１の半径及び内壁酸化膜２２の厚さを求める計算を停止し、かつ酸素析出核の半径を求める。酸素析出核は酸素と空孔を消費しながら成長し、空孔はＳｉとその酸化物（ＳｉＯｘ）との間の歪みを調整する。また本発明のモデルでは、酸素析出核の成長速度は酸素の拡散と空孔の拡散に限定すると仮定している。γＤ_OＣ_O＜Ｄ_v（Ｃ_v−Ｃ_ve）であるならば、酸素析出核の成長速度は酸素の拡散に限定され、次の式（２９）により求められる。
【数２３】

【００５３】
またγＤ_OＣ_O＞Ｄ_v（Ｃ_v−Ｃ_ve）であるならば、クラスターの成長時に、Ｓｉとその酸化物（ＳｉＯｘ）との間の歪みを無くすために空孔が消費されるので、クラスターの成長速度は次の式（３０）で表される。
【数２４】

【００５４】
上記第８ステップから第１６ステップをシリコン単結晶１４の冷却が完了するまで繰返す（第１７ステップ）。なお、式（９）〜式（３０）は互いに関連させてコンピュータにより解く。
【００５５】
上述のように、シリコン融液１２の対流を考慮してシリコン融液１２から成長するシリコン単結晶１４内の温度分布を求めた後に、シリコン融液１２から切離されたシリコン単結晶１４の冷却過程を考慮し、シリコン単結晶１４の徐冷及び急冷の効果を結果に反映して解析することにより、シリコン単結晶１４内の酸素析出核の密度分布及びサイズ分布を正確に予測できる。即ち、シリコン単結晶１４がシリコン融液１２から切離された後のシリコン単結晶１４の引上げ速度を考慮して、ボイド２１及び内壁酸化膜２２からなるボイド欠陥の密度分布とサイズ分布をコンピュータを用いて求め、酸素析出核の単位時間当りの発生量Ｊ_cとこの臨界半径ｒ_crとをコンピュータを用いて求め、更に臨界半径ｒ_cr以上の半径ｒ_cを有する酸素析出核の単位時間当りの発生量Ｊ_cがゼロを越えたときに、ボイド半径ｒ_c及び内壁酸化膜厚さｄを求める計算を停止し、かつ臨界半径ｒ_cr以上の半径ｒ_cを有する酸素析出核の半径ｒ_cをコンピュータを用いて求める。この結果、臨界半径ｒ_cr以上の半径ｒ_cを有する酸素析出核の単位時間当りの発生量Ｊ_cを求めることにより、シリコン単結晶１４内の酸素析出核の密度分布を正確に予測できる。また臨界半径ｒ_cr以上の半径ｒ_cを有する酸素析出核の半径ｒ_cを求めることにより、シリコン結晶１４内の酸素析出核のサイズ分布を正確に予測できる。
なお、この実施の形態では、単結晶としてシリコン単結晶を挙げたが、ＧａＡｓ単結晶，ＩｎＰ単結晶，ＺｎＳ単結晶若しくはＺｎＳｅ単結晶でもよい。
【００５６】
【実施例】
次に本発明の実施例を比較例とともに詳しく説明する。
＜実施例１＞
図５及び図６に示すように、石英るつぼ１３に貯留されたシリコン融液１２から直径６インチのシリコン単結晶１４を引上げる場合の、シリコン単結晶１４内のボイド２１及び内壁酸化膜２２からなるボイド欠陥の密度分布及びサイズ分布を、図１〜図４のフローチャートに基づくシミュレーション方法により求めた。
【００５７】
即ち、シリコン単結晶引上げ機１１のホットゾーンをメッシュ構造でモデル化した。ここで、シリコン融液１２のシリコン単結晶１４直下のシリコン単結晶１４の径方向のメッシュを０．７５ｍｍに設定し、シリコン融液１２のシリコン単結晶１４直下以外のシリコン単結晶１４の径方向のメッシュを１〜５ｍｍに設定した。またシリコン融液１２のシリコン単結晶１４の長手方向のメッシュを０．２５〜５ｍｍに設定し、乱流モデル式の乱流パラメータＣとして０．４５を用いた。更にシリコン単結晶１４の引上げ開始時ｔ₀から冷却完了時ｔ₁までの引上げ長及び引上げ高さの段階的な変更を５０ｍｍずつとした。このような条件下で、シリコン融液１２の対流を考慮してシリコン単結晶１４内の温度分布を求めた。
【００５８】
次にシリコン単結晶１４の引上げ開始時ｔ₀から冷却完了時ｔ₁までの時間を所定の間隔Δｔ秒（微小な時間の間隔）毎に区切り、この時間間隔Δｔ秒毎に上記シリコン単結晶１４のメッシュの座標及び温度のデータから、シリコン単結晶１４の引上げ長及び引上げ高さとシリコン単結晶内１４の温度分布を求め、更にボイド２１の形成開始温度を求めてボイド２１の密度分布及びサイズ分布を求めた。即ち、シリコン単結晶１４をシリコン融液１２から切離した後のシリコン単結晶１４の徐冷及び急冷を考慮して、シリコン単結晶１４内のボイド２１の密度分布及びサイズ分布をコンピュータを用いてそれぞれ求めた。
【００５９】
一方、ボイド２１の周囲に成長する内壁酸化膜２２は、ボイド２１が発生するとすぐに成長するものとして扱ったので、内壁酸化膜２２の形成開始温度Ｔ_pをボイド２１の形成開始温度Ｔ_vと同一とした。ここで、「すぐに」とは、「ボイドが発生した次の繰返し計算から」という意味である。次に内壁酸化膜２２の形成開始温度より低くなったときの内壁酸化膜２２の外半径Ｒ_pをコンピュータを用いて求めた。次に壁酸化膜２２の外半径Ｒ_pをボイド２１の半径から引くことにより内壁酸化膜２２の厚さｄを求めた。なお、内壁酸化膜の密度分布は上記ボイドの密度分布と同一とした。
【００６０】
次にシリコン単結晶１４内の酸素の密度分布と空孔の密度分布と温度分布に基づいて、酸素析出核の単位時間当りの発生量Ｊ_cとこの臨界半径ｒ_crとをコンピュータを用いて求めた。更に臨界半径ｒ_cr以上の半径ｒ_cを有する酸素析出核の単位時間当りの発生量Ｊ_cがゼロを越えたときに、ボイド半径ｒ_c及び内壁酸化膜厚さｄを求める計算を停止し、かつ臨界半径ｒ_cr以上の半径ｒ_cを有する酸素析出核の半径ｒ_cをコンピュータを用いて求めた。
【００６１】
シリコン単結晶１４の半径方向の位置を変えた場合（図８（ａ）のＡ部とＢ部）の、上記酸素析出核のサイズ分布及び密度分布を図８（ｂ）の実線Ａ及び実線Ｂで示した。またシリコン単結晶１４の引上げ速度を変えた場合（図８（ａ）のＡ部とＣ部）の、上記酸素析出核のサイズ分布及び密度分布を図８（ｂ）の実線Ａ及び実線Ｃで示した。
【００６２】
＜比較例１＞
実施例１と同一形状のシリコン単結晶を同一の条件で引上げた。このシリコン単結晶を比較例１とした。
【００６３】
＜比較試験及び評価＞
実施例１の方法で求めた図８（ａ）のＡ部、Ｂ部及びＣ部での酸素析出核の密度分布のうち、経験的に求めた熱処理条件に依存したあるサイズ以上の酸素析出核の総数（Ａ₁、Ｂ₁及びＣ₁）をコンピュータを用いて求めた。その結果を表１に示す。一方、比較例１のシリコン単結晶において、上記図８（ａ）のＡ部、Ｂ部及びＣ部に対応する部分の酸素析出核の密度分布を所定の熱処理後に光学顕微鏡にて測定した。その結果を表１に示す。
【００６４】
【表１】

【００６５】
表１から明らかなように、実施例１の方法で算出した酸素析出核の密度分布は比較例１の実際に測定した酸素析出核の密度分布とオーダーが一致した。この結果、酸素析出核の密度分布を実施例１の方法によりオーダーレベルで見積もり可能であることが判った。
【００６６】
【発明の効果】
以上述べたように、本発明によれば、融液の対流を考慮して融液から成長する単結晶内の温度分布をコンピュータを用いて求めるだけでなく、更に冷却過程における単結晶内の温度分布までも求めることによって、単結晶内の酸素析出核の密度分布及びサイズ分布を正確に予測できる。即ち、融液から切離された単結晶の冷却過程における単結晶の徐冷及び急冷の効果を考慮することによって、ボイド及び内壁酸化膜からなるボイド欠陥の密度分布とサイズ分布を求め、酸素析出核の単位時間当りの発生量とこの臨界半径とを求め、更に酸素析出核の単位時間当りの発生量がゼロを越えたときに、ボイド半径及び内壁酸化膜厚さを求める計算を停止し、かつ酸素析出核の半径をコンピュータを用いて求める。この結果、酸素析出核の単位時間当りの発生量を求めることにより単結晶内の酸素析出核の密度分布を正確に予測でき、酸素析出核の半径を求めることにより結晶内の酸素析出核のサイズ分布を正確に予測できる。
【図面の簡単な説明】
【図１】本発明実施形態シリコン単結晶内のボイド欠陥の密度分布及びサイズ分布のシミュレーション方法の第１段を示すフローチャート。
【図２】そのボイド欠陥の密度分布及びサイズ分布のシミュレーション方法の第２段を示すフローチャート。
【図３】そのボイド欠陥の密度分布及びサイズ分布のシミュレーション方法の第３段を示すフローチャート。
【図４】そのボイド欠陥の密度分布及びサイズ分布のシミュレーション方法の第４段を示すフローチャート。
【図５】本発明のシリコン融液をメッシュ構造としたシリコン単結晶の引上げ機の要部断面図。
【図６】そのシリコン単結晶内のボイド及び内壁酸化膜の模式図。
【図７】従来例のシリコン融液をメッシュ構造としたシリコン単結晶の引上げ機の要部断面図。
【図８】シリコン単結晶の半径方向の位置を変えた場合と、シリコン単結晶の引上げ速度を変えた場合に、コンピュータを用いて求めた酸素析出核のサイズ分布及び密度分布を示す図。
【符号の説明】
１１シリコン単結晶引上げ機
１２シリコン融液
１４シリコン単結晶
２１ボイド
２２内壁酸化膜
Ｓシリコンの三重点[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a method for computer simulation of the density distribution and size distribution of oxygen precipitation nuclei in a single crystal such as silicon pulled by the Czochralski (hereinafter referred to as CZ) method. The oxygen precipitation nuclei mentioned here are oxygen precipitation nuclei formed by consuming vacancies and oxygen in a single crystal, and this is stabilized by a low-temperature (700 to 800 ° C.) heat treatment, and further high temperature (1000 to When grown by heat treatment at 1100 ° C., it is observed as BMD (Bulk Micro-Defect) with an optical microscope or the like.
[0002]
[Prior art]
Conventionally, as a simulation method of this type, as shown in FIG. 7, the hot zone structure in the pulling machine 1 and the pulling speed of the silicon single crystal 4 when pulling the silicon single crystal 4 by the CZ method using a general heat transfer simulator. Based on the above, the internal temperature distribution of the silicon melt 2 is predicted by manipulating the thermal conductivity of the silicon melt 2, and the mesh coordinates and temperature of the silicon single crystal 4 are obtained from the internal temperature distribution, respectively. A method is known in which the density distribution of the interstitial silicon and vacancies is obtained by using a computer by solving the diffusion equation based on the diffusion coefficient and boundary conditions of the interstitial silicon and vacancies in the single crystal 4. In this simulation method, each member of the hot zone is divided into meshes and modeled. In particular, the mesh of the silicon melt 2 is set to be relatively coarse, such as about 10 mm, in order to shorten the calculation time.
[0003]
On the other hand, the operating conditions and the temperature conditions in the furnace when the silicon single crystal is grown by the CZ method are input as data to a computer, and the size and density of vacancy-induced grow-in defects introduced into the single crystal and the growth rate of the single crystal. Is determined by computer simulation, and the growth rate is selected so that the size and density of the vacancy-induced grow-in defects introduced into the single crystal become a predetermined value based on the simulation result. A method for producing a silicon single crystal that grows silicon is known (see, for example, Patent Document 1). In this silicon single crystal manufacturing method, the heat distribution in the furnace is simulated using a comprehensive heat transfer analysis program to determine the temperature gradient in the vertical direction of the single crystal center during single crystal growth.
[0004]
[Patent Document 1]
JP 2002-145696 A
[0005]
[Problems to be solved by the invention]
However, in the above conventional method of simulating the density distribution of interstitial silicon and vacancies, the actual puller does not consider the convection of the silicon melt generated, and the mesh of the silicon melt is relatively coarse. The reproducibility of the solid-liquid interface shape is poor and it is impossible to provide an accurate temperature distribution in the single crystal. Therefore, the density distribution of interstitial silicon and vacancies in the silicon single crystal is significantly different from the measured value, and oxygen precipitation nuclei that are defects in the silicon single crystal generated at a relatively low temperature during the crystal cooling process. The density distribution and size distribution were not known.
On the other hand, in the method for manufacturing a silicon single crystal disclosed in Patent Document 1, the size distribution and density distribution of the grow-in defects in the vertical direction of the silicon single crystal can be predicted. In addition, the density distribution cannot be predicted, and the glow-in defect on the inner surface of the wafer cannot be examined.
[0006]
The object of the present invention is to analyze the temperature distribution in a growing single crystal in consideration of the convection of the melt and then analyze the cooling process of the single crystal separated from the melt. An object is to provide a method for simulating the density distribution and size distribution of oxygen precipitation nuclei in a single crystal, which can accurately predict the density distribution and size distribution of oxygen precipitation nuclei in the crystal.
[0007]
[Means for Solving the Problems]
As shown in FIGS. 1 to 6, the invention according to claim 1 is a hot zone of the puller 11 from the start of pulling of the single crystal 14 from the melt 12 by the puller 11 to the completion of cooling of the single crystal 14. The first step of modeling the mesh with a mesh structure, the mesh for each member of the hot zone, and the single crystal 14 corresponding to the pulling length of the single crystal 14 together with the physical properties of each member with respect to the combined mesh A second step of inputting the pulling speed of 14 to the computer, a third step of obtaining the surface temperature distribution of each member based on the amount of heat generated by the heater and the emissivity of each member, and the surface temperature distribution and heat conduction of each member. The turbulent flow model equation 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 rate and Navier A fourth step for further obtaining the internal temperature distribution of the melt 12 in consideration of convection by connecting and solving the Stokes equation; and the solid-liquid interface shape of the single crystal 14 and the melt 12 includes the triple point S of the single crystal. The fifth step obtained according to the isotherm and the third to fifth steps are repeatedly calculated 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 determine the mesh coordinates of the single crystal 14. The sixth step of obtaining the temperature and inputting these data to the computer, respectively, and the temperature in the puller 11 are repeated from the first step to the sixth step by changing the pulling length and pulling height of the single crystal 14 stepwise. The distribution is calculated to obtain the mesh coordinates and temperature of the single crystal 14 and input these data to the computer, respectively, and the single crystal 14 is drawn from the melt 12. The time from the start of bending to the completion of cooling of the single crystal 14 is divided at predetermined intervals, and the single crystal 14 is obtained from the mesh coordinates and temperature data of the single crystal 14 obtained in the seventh step at each divided time interval. The eighth step for obtaining the pulling length and pulling height and the temperature distribution in the single crystal 14, and solving the diffusion equation based on the diffusion coefficients and boundary conditions of vacancies and interstitial atoms in the single crystal 14 A ninth step for obtaining the density distribution of vacancies and interstitial atoms after the elapse of the time interval, a tenth step for obtaining the formation start temperature of the void 21 based on the density distribution of the vacancies, and the single crystal 14 respectively. An eleventh step of determining the density of the void 21 when the temperature at the lattice point of the mesh gradually decreases to the formation start temperature of the void 21, and each mesh in the single crystal 14 A twelfth step for determining the radius of the void 21 when the temperature at the lattice point is lower than the formation start temperature of the void 21, a thirteenth step for determining the formation start temperature of the inner wall oxide film 22 grown around the void 21, The radius of the void 21 and the thickness of the inner wall oxide film 22 when the temperature at the lattice point of each mesh in the crystal 14 is gradually lowered to be lower than the formation start temperature of the inner wall oxide film 22 are obtained in relation to each other. A fifteenth step, a fifteenth step for determining the amount of oxygen precipitation nuclei generated per unit time and a critical radius based on the density distribution of oxygen in the single crystal 14, the density distribution of vacancies, and the temperature distribution; When the generation amount per unit time exceeds zero, the calculation for determining the radius of the void 21 and the thickness of the inner wall oxide film 22 is stopped and the sixteenth step for determining the radius of the oxygen precipitation nucleus is performed. And a simulation of the density distribution and the size distribution of the oxygen precipitation nuclei in the single crystal using a computer including the eighteenth step to the sixteenth step until the cooling of the single crystal 14 is completed. .
[0008]
In the simulation method according to the first aspect, not only the temperature distribution in the single crystal 14 grown from the melt 12 is obtained in consideration of the convection of the melt 12, but also the temperature in the single crystal 14 in the cooling process. By determining even the distribution, the density distribution and size distribution of the oxygen precipitation nuclei in the single crystal 14 can be accurately predicted. That is, by analyzing the effect of slow cooling and rapid cooling of the single crystal 14 in the cooling process of the single crystal 14 separated from the melt 12, the density distribution of void defects composed of the void 21 and the inner wall oxide film 22 is analyzed. And the size distribution, the amount of oxygen precipitation nuclei generated per unit time and the critical radius are obtained, and when the amount of oxygen precipitation nuclei generated per unit time exceeds zero, the void radius and the inner wall oxide film thickness Is stopped, and the radius of the oxygen precipitation nuclei is obtained. Thereby, the density distribution and size distribution of oxygen precipitation nuclei in the single crystal 14 can be accurately predicted.
Here, the critical radius of the oxygen precipitation nucleus means the minimum radius of the nucleus generated during crystal growth, and this radius depends on the supersaturation degree of the vacancy density and the supersaturation degree of the oxygen density.
[0009]
DETAILED DESCRIPTION OF THE INVENTION
Next, embodiments of the present invention will be described with reference to the drawings.
As shown in FIG. 5, a quartz crucible 13 for storing the silicon melt 12 is provided in the chamber of the silicon single crystal puller 11. Although not shown, the quartz crucible 13 is connected to a crucible driving means via a graphite susceptor and a support shaft, and the crucible driving means is configured to rotate and raise and lower the quartz crucible 13. Further, the outer peripheral surface of the quartz crucible 13 is surrounded by a heater (not shown) at a predetermined interval from the quartz crucible 13, and the heater is surrounded by a heat insulating cylinder (not shown). The heater heats and melts the high-purity silicon polycrystal charged in the quartz crucible 13 to form the silicon melt 12. A cylindrical casing (not shown) is connected to the upper end of the chamber, and this casing is provided with a pulling means. The pulling means is configured to pull the silicon single crystal 14 while rotating it.
[0010]
A method of simulating the density distribution and size distribution of oxygen precipitation nuclei in the silicon single crystal 14 in the silicon single crystal puller 11 configured as described above will be described with reference to FIGS.
First, as a first step, a silicon single crystal 14 has a predetermined length L.₁Each member of the hot zone of the silicon single crystal pulling machine 11 in a state of being pulled up (for example, 100 mm), that is, the chamber, the quartz crucible 13, the silicon melt 12, the silicon single crystal 14, the graphite susceptor, the heat insulating cylinder, etc. are divided into meshes. Model. Specifically, coordinate data of mesh points of each member in the hot zone is input to the computer. At this time, among the mesh of the silicon melt 12, a mesh in the radial direction of the silicon single crystal 14 and a part or all of the mesh immediately below the silicon single crystal 14 of the silicon melt 12 (hereinafter referred to as a radial mesh). The thickness is set to 0.01 to 5.00 mm, preferably 0.25 to 1.00 mm. Of the mesh of the silicon melt 12, a mesh in the longitudinal direction of the silicon single crystal 14 and a part or all of the silicon melt 12 (hereinafter referred to as a longitudinal mesh) is 0.01 to 5.00 mm. The thickness is preferably set to 0.1 to 0.5 mm.
[0011]
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 silicon melt 12 near the outer peripheral edge of the silicon single crystal 14 out of the silicon melt 12 immediately below the silicon single crystal 14 is used. It is preferable to limit to the above range. When a part of the longitudinal mesh is limited to the range of 0.01 to 5.00, the vicinity of the liquid surface and the bottom of the silicon melt 12 is limited to the above range. Is preferred.
[0012]
As a second 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 with respect to 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 single crystal 14, the pulling speed of the silicon single crystal 14 corresponding to the pulling length, and the turbulent flow parameter C of the turbulent flow model equation (1) described later are input to the computer.
[0013]
As a third step, the surface temperature distribution of each member in the hot zone is obtained 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 inputted to the computer, and the surface temperature distribution of each member is obtained from the radiation rate of each member using the computer. Next, as a fourth step, the internal temperature distribution of each member is obtained 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, the xyz orthogonal coordinate system is used, but in the actual calculation, a cylindrical coordinate system is used.
[Expression 1]

[0014]
In the above formula (1), ρ is the density of each member, c is the specific heat of each member, T is the absolute temperature at each mesh point of each member, and t is time. Λ_x, Λ_yAnd λ_zIs the x, y and z direction components of the thermal conductivity of each member, and q is the amount of heat generated by the heater.
[0015]
On the other hand, for the silicon melt 12, after obtaining the internal temperature distribution of the silicon melt 12 by the above heat conduction equation (1), the silicon melt 12 is turbulent based on the internal temperature distribution of the silicon melt 12. By connecting the turbulent flow model equation (2) and the Navier-Stokes equations (3) to (5) obtained by assuming that there is, the internal flow velocity distribution of the silicon melt 12 is obtained using a computer.
[Expression 2]

[0016]
In the above formula (2), κ_tIs the turbulent thermal conductivity of the silicon melt 12, c is the specific heat of the silicon melt 12, Pr_tIs the Prandtl number, ρ is the density of the silicon melt 12, C is the turbulent parameter, d is the distance from the quartz crucible 13 wall storing the silicon melt 12, and k is the silicon melt 12. The sum of squares of the fluctuation component with respect to the average flow velocity.
[Equation 3]

[0017]
In the above formulas (3) to (5), u, v, and w are x, y, and z direction components of the flow velocity at each mesh point of the silicon melt 12, and v_lIs the molecular kinematic viscosity coefficient (physical property value) of the silicon melt 12 and v_tIs the kinematic viscosity coefficient due to the effect of turbulent flow of the silicon melt 12, and F_x, F_yAnd F_zAre the x, y and z direction components of the body force acting on the silicon melt 12.
[0018]
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 silicon melt 12 is a fluid that is incompressible and has a constant viscosity.
[0019]
By solving the thermal energy equation (6) based on the obtained internal flow velocity distribution of the silicon melt 12, the internal temperature distribution of the silicon melt 12 considering the convection of the silicon melt 12 is further obtained using a computer. .
[Expression 4]

[0020]
In the above equation (6), u, v and w are x, y and z direction components of the flow velocity at each mesh point of the silicon melt 12, and T is an absolute temperature at each mesh point of the silicon melt 12. Ρ is the density of the silicon melt 12, c is the specific heat of the silicon melt 12, κ_lIs the molecular thermal conductivity (physical property value) and κ_tIs the turbulent thermal conductivity calculated using equation (1).
[0021]
Next, as a fifth step, a silicon triple point S (solid-liquid-gas triple point (tri-junction)) in which the solid-liquid interface shape of the silicon single crystal 14 and the silicon melt 12 is indicated by a point S in FIG. 5 is included. Use a computer to match the isotherm. As a sixth step, the amount of heat generated by the heater input to the computer is changed (increase gradually), and the third to fifth steps are repeated until the triple point reaches the melting point of the silicon single crystal 14, and then the puller 11 The temperature distribution is calculated to determine the coordinates and temperature of the silicon single crystal mesh, and these data are stored in the computer.
[0022]
Next, as a seventh step, the pulling length L of the silicon single crystal 14₁Δ (for example, 50 mm) is added to the above and the above first to sixth steps are repeated, and then the temperature distribution in the puller 11 is calculated to obtain the coordinates and temperature of the mesh of the silicon single crystal 14, and these data Is stored in the computer. In this seventh step, the pulling length L of the silicon single crystal 14 is₁Is length L₂(L₂Is the length of the silicon single crystal 14 when separated from the silicon melt 12 (the crystal length when the growth is completed). ) And the silicon single crystal 14 is separated from the silicon melt 12, and then the silicon single crystal 14 is further pulled up to its height H₁(H₁Is the distance from the straight cylinder start part of the silicon single crystal 14 to the liquid surface of the silicon melt 12 (FIG. 5). ) Is H₂(H₂Is the distance from the straight cylinder start part of the silicon single crystal 14 to the liquid surface of the silicon melt 12 when cooling is completed. ) Until the cooling of the silicon single crystal 14 is completed. After the silicon single crystal 14 is separated from the silicon melt 12, the pulling height H of the silicon single crystal 14 is increased.₁Δ (for example, 50 mm) is added to the above, and the first to sixth steps are repeated in the same manner as described above.
[0023]
Lifting height H of silicon single crystal 14₁Is H₂When reached, the process proceeds to the eighth step. In the eighth step, when the silicon single crystal 14 is grown from the silicon melt 12 and started to be pulled, t₀Then, the silicon single crystal 14 is separated from the silicon melt 12, and the silicon single crystal 14 is further pulled up.₁Is divided every predetermined interval Δt seconds (minute time interval). At this time, not only the diffusion coefficients and boundary conditions of interstitial silicon and vacancies in the silicon single crystal 14 but also the equations for obtaining the density distribution and size distribution of the void 21 and the inner wall oxide film 22 (FIG. 6) described later. Each constant to be entered is input to the computer. The pulling length L of the silicon single crystal 14 is obtained from the mesh coordinates and temperature data of the silicon single crystal 14 obtained in the seventh step at every divided time interval Δt seconds.₁And pulling height H₁And the temperature distribution in the silicon single crystal 14 are obtained.
[0024]
That is, in the first to seventh steps, the mesh coordinates and temperature of the silicon single crystal 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 t₀The pulling length L of the silicon single crystal 14 after Δt seconds from₁And pulling height H₁And the temperature distribution in the silicon single crystal 14 is calculated. Next, the density distribution of the vacancies and interstitial silicon after lapse of Δt seconds is obtained by solving the diffusion equation based on the diffusion coefficients and boundary conditions of the vacancies and interstitial silicon in the silicon single crystal 14 (9th step). ).
[0025]
Specifically, the hole density C_vIs calculated by the following formula (7), and the density C of interstitial silicon is_iIs calculated by the following equation (8). In the equations (7) and (8), the density C_vAnd density C_iIn 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.
[Equation 5]

[0026]
In the above formulas (7) and (8), K₁And K₂Is a constant, E_iAnd E_vAre the energy of formation of interstitial silicon and vacancies, respectively, and C_veAnd C_ieIs the equilibrium density of vacancies and the equilibrium density of interstitial silicon. K_BMeans Boltzmann constant and T means absolute temperature.
[0027]
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.
[Formula 6]

[0028]
In the above formulas (9) and (10), Θ (x) is a heaviside function. That is, Θ (x) = 0 when x <0, and Θ (x) = 1 when x> 0. T_pIs the formation start temperature of the inner wall oxide film 22 (silicon oxide film (SiOx film) formed on the surface of the void 21), and T_vIs the formation start temperature of the void 21 (FIG. 6). D_OIs the oxygen diffusion coefficient, C_OIs the oxygen density at the oxygen precipitation nucleation temperature. J_cIs the amount of oxygen precipitation nuclei generated per unit time with a critical radius, r_cIs the radius of the oxygen precipitation nuclei. The first term on the right side of each of the equations (9) and (10) is Fick's diffusion formula, and D in the first term on the right side is D_vAnd D_iIs a diffusion coefficient of vacancies and interstitial silicon expressed by the following equations (11) and (12).
[Expression 7]

[0029]
In the above formula (11) and formula (12), ΔE_vAnd △ E_iAre the activation energies of the vacancies and interstitial silicon, respectively, d_vAnd d_iAre constants. Further, E in the second term on the right side of each of the equations (9) and (10)_v ^tAnd E_i ^tIs the activation energy of vacancies and interstitial silicon due to thermal diffusion, k_BIs the Boltzmann constant. The third term k on the right side of each of the equations (9) and (10)_ivIs the recombination constant of vacancies and interstitial silicon pairs. N in the fourth term on the right side of each of the equations (9) and (10)_vIs the density of void 21 and r_vIs the radius of the void 21, and N in the fifth term on the right side of the equation (9)_pIs the density of the inner wall oxide film 22 and R_pIs the outer radius of the inner wall oxide film 22.
[0030]
The above equation (9) is satisfied because the vacancy flux for vacancy deposition is sufficiently large, and the unit mass of the Si matrix constituting the silicon single crystal 14 and the SiOx constituting the inner wall oxide film per unit mass. When the volume difference can be filled, that is, γD_OC_O≦ D_v(C_v-C_ve). In other cases, that is, γD_OC_O> D_v(C_v-C_ve), The following equation (13) is established. Here, γ is the amount of vacancies per oxygen atom necessary for the oxide to grow without distortion, and is a number depending on the phase of the oxygen precipitation nuclei. Oxygen precipitation nuclei are SiO₂When γ is approximately γ≈0.5, and when the oxygen precipitation nucleus is SiO, γ is approximately 0.65. D_OIs the diffusion coefficient of oxygen, C_OIs the oxygen density at the oxygen precipitation nucleation temperature.
[Equation 8]

[0031]
In the above equation (13), Θ (x) is a heaviside function. That is, Θ (x) = 0 when x <0, and Θ (x) = 1 when x> 0. T_pIs the temperature at which the inner wall oxide film 22 is formed, and T_vIs the formation start temperature of the void 21. J_cIs the amount of oxygen precipitation nuclei generated per unit time with a critical radius, r_cIs the radius of the oxygen precipitation nuclei. The first term on the right side of equation (13) is Fick's diffusion formula, and D in the first term on the right side_vIs a diffusion coefficient of holes represented by the formula (11).
[0032]
Next, as a tenth step, the density C of the pores obtained by solving the above diffusion equation_vBased on the distribution, formation start temperature T of void 21_vIs obtained from the following equation (14).
[Equation 9]

[0033]
In the above formula (14), C_vIs the vacancy density in the silicon single crystal 14 and C_v0Is the melting point T of the silicon single crystal 14_mIs the vacancy equilibrium density at E_vIs the vacancy formation energy. Also σ_vIs the interfacial energy at the crystal plane (111) of the silicon single crystal 14, ρ is the density of the silicon single crystal 14, and k_BIs the Boltzmann constant.
[0034]
As an eleventh step, the temperature at the lattice point of each mesh in the silicon single crystal 14 gradually decreases, and the formation start temperature T of the void 21 is increased._vThe density N of the void 21 using the following approximate expression (15)_vAsk for.
[Expression 10]

[0035]
In the above formula (15), n_vIs the density of the void 21 and (dT / dt) is the cooling rate of the silicon single crystal 14. Although the density of the void 21 depends on the vacancy supersaturation degree, it may be considered constant in a narrow temperature region such as the nucleation temperature in the silicon single crystal 14 to be calculated. D_vIs the diffusion coefficient of void 21 and k_BIs the Boltzmann constant and C_vIs the vacancy density in the silicon single crystal 14.
[0036]
As a twelfth step, the temperature at the lattice point of each mesh in the silicon single crystal 14 is the formation start temperature T of the void 21._vLower void radius r_vIs obtained from the following equation (16). Here, the growth rate of the void 21 is a diffusion-controlled rate that depends on the diffusion rate of the holes. Although the shape of the void 21 is actually an octahedron, it is treated as a sphere for efficiency of calculation here.
## EQU11 ##

[0037]
In the above equation (16), t1 is the temperature at the lattice point of the mesh of the silicon single crystal 14 and the formation start temperature T of the void 21_vIt is the time when it drops to. R_vcrIs the critical radius of the void 21 and this critical radius r_vcrIs the value at time t1. In addition C_veAnd C_ieIs the equilibrium density of vacancies and interstitial silicon in each state.
[0038]
As a thirteenth step, the formation start temperature T of the inner wall oxide film 22 growing around the void 21_pAsk for. Formation start temperature T of inner wall oxide film 22_pIs the formation start temperature T of the void 21_vIf it is smaller, oxygen and vacancies combine to grow an inner wall oxide film 22 on the surface of the void 21 (FIG. 6). In the model of the present invention, the inner wall oxide film 22 is treated as growing immediately after the void 21 is generated. Here, “immediately” means “from the next repeated calculation in which a void has occurred”.
[0039]
As a fourteenth step, the temperature at the lattice point of each mesh in the silicon single crystal 14 gradually decreases, and the formation start temperature T of the inner wall oxide film 22 is decreased._pThe radius r of the void 21 when it becomes lower_vAnd the thickness d of the inner wall oxide film 22 are obtained in relation to each other. Generally, the inner wall oxide film 22 is formed by the diffusion of oxygen and the reaction with silicon. Flux J in which oxygen flows from outside the inner wall oxide film 22 to the outer peripheral surface of the oxide film_OR denotes the outer radius of the inner wall oxide film 22 with the center of the void 21 as the origin._pIs obtained from the following equation (17).
[Expression 12]

[0040]
In the above formula (17), D_OIs the diffusion coefficient of oxygen, C_OIs the oxygen density and C_OeifIs the radius R with the center of the spherical void 21 as the origin_pIs the oxygen equilibrium density in g_OIs the reaction rate of oxygen atoms at the interface between SiOx and Si, that is, the outer peripheral surface of the inner wall oxide film 22. In addition, the voids are formed from outside the inner wall oxide film 22 to the outer peripheral surface of the oxide film (radius R_pFlux flowing into the spherical surface_vIs obtained from the following equation (18).
[Formula 13]

[0041]
In the above formula (18), g_vIs the reaction rate of vacancies on the outer peripheral surface of the inner wall oxide film 22, and C_veifIs the equilibrium density of holes on the outer peripheral surface of the inner wall oxide film 22. As the inner wall oxide film 22 of silicon grows, vacancies are consumed, and the vacancy consumption rate is γ for one oxygen atom. Therefore, the condition for the inner wall oxide film 22 to grow without distortion is the oxygen flux J_OIn contrast, the hole flux is γJ_OSo the flux of holes γJ_OIs obtained from the following equation (19).
[Expression 14]

[0042]
Therefore, (J_v-ΓJ_O) Is a flux of holes absorbed by the inner wall oxide film 22. Here, when the flux from outside the inner wall oxide film 22 of the pores to this oxide film increases ((J_v-ΓJ_O)> 0), or when the flux from the void 21 to the oxide film increases with time ((J_v-ΓJ_O) <0) where the rate of change of the flux with respect to time is α, this α is considered to depend on the thickness of the inner wall oxide film 22. Therefore, the thickness d of the inner wall oxide film 22, that is, (R_p-R_v) Is d₀If it is less than that, it is considered that the inner wall oxide film 22 partially penetrates the void from both sides of the oxide film and the void 21. In other words, d <d₀If α ≠ 0, vacancies penetrate the inner wall oxide film 22 at a certain rate and are used for the growth of the void 21. However, the interstitial silicon is not used for the growth of the void 21 at the same time as the oxygen precipitation layer is nucleated.
[0043]
Further, when the amount of flux from outside the inner wall oxide film 22 to the oxide film is small, that is, (J_v-ΓJ_O) <0, the inner wall oxide film 22 without distortion cannot be grown, so that the voids flow into the oxide film from the void 21 and the outer diameter of the void 21 is reduced. As a result, d ≧ d₀When α = 0, the inner wall oxide film 22 does not pass through the holes at all. And n_vWhere dn is the number of vacancies in void 21_v/ dt is the inflow of holes, dn_v/ dt = α (J_v-ΓJ_O), The following equation (20) holds. In formula (20), r_vIs the radius of the void 21 obtained by the equation (16), and R_pIs the outer radius of the inner wall oxide film 22. D_OIs the diffusion coefficient of oxygen, C_OIs the oxygen density. In addition C_OeifIs the radius R with the center of the spherical void 21 as the origin_pIs the oxygen equilibrium density in g_OIs the reaction rate of oxygen atoms on the outer peripheral surface of the inner wall oxide film 22.
[Expression 15]

[0044]
At the same time, the number of SiOx molecules constituting the inner wall oxide film 22 is the flux J in which the oxygen in the formula (20) flows from the inner wall oxide film 22 to the outer peripheral surface of this oxide film._OTherefore, the following equation (21) is established. In equation (21), ρ_p= X / Ω_SiOxIt is defined as This ρ_pIn the definition formula, Ω_SiOxIs the molecular volume of SiOx, the unit of which is molecular weight / density. X corresponds to x of SiOx.
[Expression 16]

[0045]
Substituting the above equation (21) into the above equation (20) yields the following equation (22). In Equation (22), η is Ω_SiOx/ Ω_SiIt is. In this definition of η, Ω_SiIs 1 / ρ, and ρ is the density of the silicon single crystal 14.
[Expression 17]

[0046]
The following can be understood from the above equations (20) and (22). First, d is d₀T_ThreeThen, this time t_ThreeIn this case, α becomes zero, so that the inner wall oxide film 22 does not penetrate. The outer diameter of the void 21 is r_v(t_Three), The growth of the inner wall oxide film 22 is determined by the consumption of vacancies or oxygen flowing from the matrix of the silicon single crystal 14. That is, (J_v-ΓJ_O)> 0, equation (22) becomes the following equation (23), and (J_v-ΓJ_O) <0, Equation (22) becomes the following Equation (24).
[Expression 18]

[0047]
In the above formula (23) and formula (24), Ω_SiOxIs the molecular volume of SiOx. In the equations (23) and (24), the vacancy equilibrium density on the outer peripheral surface of the inner wall oxide film is R_pIs unrelated to C_veif= C_veAnd C_Oeif= C_OeIs assumed. From equation (23) or equation (24), outer radius R of inner wall oxide film 22_pThe outer radius R of the inner wall oxide film 22_pTo radius 21 of void 21_vTo obtain the thickness d of the inner wall oxide film 22.
[0048]
Next, as a fifteenth step, based on the oxygen density distribution, the vacancy density distribution, and the temperature distribution in the silicon single crystal 14, the amount of oxygen precipitation nuclei generated per unit time J_cAnd the critical radius r_crAnd ask. Specifically, oxygen precipitation nucleus growth (SiOx phase) is performed by oxygen density C_OAnd pore density C_vIn the classical theory, the amount of oxygen precipitation nuclei generated per unit time in the steady state is J_cIs obtained from the following equation (25).
[Equation 19]

[0049]
In the above equation (25), σ_cIs the interface energy between SiOx and Si at the oxygen precipitation nucleation temperature, D_OIs the oxygen diffusion coefficient. Also C_OIs the oxygen density at the oxygen precipitation nucleation temperature and k_BIs the Boltzmann constant. F^*Is a growth barrier for oxygen precipitation nuclei, and assuming that the oxygen precipitation nuclei are spheres, F^*Can be obtained by the following equation (26).
[Expression 20]

[0050]
In the above equation (26), ρ_cIs the oxygen number density of the oxygen precipitation nuclei made of SiOx, and f_cIs the driving force (force per oxygen atom) necessary for the formation of oxygen precipitation nuclei, and this f_cIs obtained from the following equation (27).
[Expression 21]

[0051]
In the above formula (27), C_OIs the oxygen density at the oxygen precipitation nucleation temperature, C_OeIs the oxygen equilibrium density. Γ is the amount of vacancies per one oxygen atom necessary for the oxide to grow without distortion, and is a number depending on the phase of the oxygen precipitation nuclei. Oxygen precipitation nuclei are SiO₂When γ is approximately γ≈0.5, and when the oxygen precipitation nucleus is SiO, γ is approximately 0.65. In addition C_vIs the density of vacancies and C_veIs the equilibrium density of vacancies. Note that σ in Equation (25) and Equation (26)_c(Interfacial energy between SiOx and Si at the oxygen precipitation nucleation temperature) is an important parameter governing nucleation, and this parameter σ_cDepends on nucleation, and generally varies depending on temperature and radius of oxygen precipitation nuclei. On the other hand, the critical radius r of the oxygen precipitation nucleus_crIs obtained by the following equation (28).
[Expression 22]

[0052]
Next, as a sixteenth step, when the generation amount of oxygen precipitation nuclei per unit time exceeds zero, the calculation for determining the radius of the void 21 and the thickness of the inner wall oxide film 22 is stopped and the radius of the oxygen precipitation nuclei is obtained. . The oxygen precipitation nuclei grow while consuming oxygen and vacancies, and the vacancies adjust the strain between Si and its oxide (SiOx). In the model of the present invention, it is assumed that the growth rate of oxygen precipitation nuclei is limited to oxygen diffusion and vacancy diffusion. γD_OC_O<D_v(C_v-C_ve), The growth rate of the oxygen precipitation nuclei is limited to the diffusion of oxygen, and is obtained by the following equation (29).
[Expression 23]

[0053]
ΓD_OC_O> D_v(C_v-C_ve), Vacancies are consumed during cluster growth to eliminate strain between Si and its oxide (SiOx), and the cluster growth rate is expressed by the following equation (30). The
[Expression 24]

[0054]
The eighth to sixteenth steps are repeated until the cooling of the silicon single crystal 14 is completed (17th step). Equations (9) to (30) are related to each other and solved by a computer.
[0055]
As described above, after determining the temperature distribution in the silicon single crystal 14 grown from the silicon melt 12 in consideration of the convection of the silicon melt 12, the silicon single crystal 14 separated from the silicon melt 12 is cooled. The density distribution and size distribution of oxygen precipitation nuclei in the silicon single crystal 14 can be accurately predicted by considering the process and analyzing the effect of slow cooling and rapid cooling of the silicon single crystal 14 in the results. That is, in consideration of the pulling speed of the silicon single crystal 14 after the silicon single crystal 14 is separated from the silicon melt 12, the density distribution and size distribution of the void defects composed of the void 21 and the inner wall oxide film 22 are calculated by a computer. The amount of oxygen precipitation nuclei generated per unit time J_cAnd this critical radius r_crUsing a computer, and the critical radius r_crRadius r above_cGeneration amount of oxygen precipitation nuclei with unit J per unit time_cVoid radius r when is over zero_cAnd the calculation of the inner wall oxide film thickness d is stopped, and the critical radius r_crRadius r above_cRadius r of oxygen precipitation nuclei having_cIs determined using a computer. As a result, the critical radius r_crRadius r above_cGeneration amount of oxygen precipitation nuclei with unit J per unit time_cThus, the density distribution of oxygen precipitation nuclei in the silicon single crystal 14 can be accurately predicted. The critical radius r_crRadius r above_cRadius r of oxygen precipitation nuclei having_cTherefore, the size distribution of the oxygen precipitation nuclei in the silicon crystal 14 can be accurately predicted.
In this 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.
[0056]
【Example】
Next, examples of the present invention will be described in detail together with comparative examples.
<Example 1>
As shown in FIGS. 5 and 6, when the silicon single crystal 14 having a diameter of 6 inches is pulled from the silicon melt 12 stored in the quartz crucible 13, the void 21 and the inner wall oxide film 22 in the silicon single crystal 14 are removed. The density distribution and size distribution of the void defects to be obtained were obtained by a simulation method based on the flowcharts of FIGS.
[0057]
That is, the hot zone of the silicon single crystal puller 11 was modeled with a mesh structure. Here, the mesh in the radial direction of the silicon single crystal 14 immediately below the silicon single crystal 14 in the silicon melt 12 is set to 0.75 mm, and the radial direction of the silicon single crystal 14 other than directly below the silicon single crystal 14 in the silicon melt 12 is set. The mesh was set to 1-5 mm. Further, the longitudinal mesh of the silicon single crystal 14 of the silicon melt 12 was set to 0.25 to 5 mm, and 0.45 was used as the turbulent flow parameter C of the turbulent flow model equation. Further, when the pulling of the silicon single crystal 14 starts₀To t when cooling is complete₁The stepwise change in the pulling length and the pulling height was set to 50 mm. Under such conditions, the temperature distribution in the silicon single crystal 14 was determined in consideration of the convection of the silicon melt 12.
[0058]
Next, when the pulling of the silicon single crystal 14 starts₀To t when cooling is complete₁Is divided at predetermined intervals Δt seconds (minute time intervals), and at each time interval Δt seconds, the pulling length of the silicon single crystal 14 and the data of the mesh coordinates and temperature of the silicon single crystal 14 are determined. The pulling height and the temperature distribution in the silicon single crystal 14 were determined, and the formation start temperature of the void 21 was further determined to determine the density distribution and size distribution of the void 21. That is, considering the slow cooling and rapid cooling of the silicon single crystal 14 after the silicon single crystal 14 is separated from the silicon melt 12, the density distribution and size distribution of the voids 21 in the silicon single crystal 14 are respectively determined using a computer. Asked.
[0059]
On the other hand, since the inner wall oxide film 22 growing around the void 21 is treated as growing immediately after the void 21 is generated, the formation start temperature T of the inner wall oxide film 22 is increased._pThe formation start temperature T of the void 21_vIt was the same. Here, “immediately” means “from the next repeated calculation in which a void has occurred”. Next, the outer radius R of the inner wall oxide film 22 when it becomes lower than the formation start temperature of the inner wall oxide film 22._pWas determined using a computer. Next, the outer radius R of the wall oxide film 22_pIs subtracted from the radius of the void 21 to obtain the thickness d of the inner wall oxide film 22. The density distribution of the inner wall oxide film was the same as the density distribution of the voids.
[0060]
Next, based on the density distribution of oxygen, the density distribution of vacancies, and the temperature distribution in the silicon single crystal 14, the amount of oxygen precipitation nuclei generated per unit time J_cAnd this critical radius r_crWas obtained using a computer. Furthermore, critical radius r_crRadius r above_cGeneration amount of oxygen precipitation nuclei with unit J per unit time_cVoid radius r when is over zero_cAnd the calculation of the inner wall oxide film thickness d is stopped, and the critical radius r_crRadius r above_cRadius r of oxygen precipitation nuclei having_cWas determined using a computer.
[0061]
When the position of the silicon single crystal 14 in the radial direction is changed (part A and part B in FIG. 8A), the size distribution and density distribution of the oxygen precipitation nuclei are shown by the solid line A and the solid line B in FIG. It showed in. Further, when the pulling speed of the silicon single crystal 14 is changed (part A and part C in FIG. 8A), the size distribution and density distribution of the oxygen precipitation nuclei are shown by a solid line A and a solid line C in FIG. Indicated.
[0062]
<Comparative Example 1>
A silicon single crystal having the same shape as in Example 1 was pulled up under the same conditions. This silicon single crystal was designated as Comparative Example 1.
[0063]
<Comparison test and evaluation>
Of the density distribution of oxygen precipitation nuclei in part A, B and C in FIG. 8A obtained by the method of Example 1, oxygen precipitation nuclei of a certain size or more depending on the heat treatment conditions obtained empirically. Total number (A₁, B₁And C₁) Using a computer. The results are shown in Table 1. On the other hand, in the silicon single crystal of Comparative Example 1, the density distribution of oxygen precipitation nuclei in the portions corresponding to the A part, the B part, and the C part in FIG. 8A was measured with an optical microscope after a predetermined heat treatment. The results are shown in Table 1.
[0064]
[Table 1]

[0065]
As is apparent from Table 1, the density distribution of oxygen precipitation nuclei calculated by the method of Example 1 was in the same order as the density distribution of oxygen precipitation nuclei actually measured in Comparative Example 1. As a result, it was found that the density distribution of oxygen precipitation nuclei can be estimated at the order level by the method of Example 1.
[0066]
【The invention's effect】
As described above, according to the present invention, not only can the temperature distribution in the single crystal grown from the melt be taken into account by taking into account the convection of the melt, but also the temperature in the single crystal during the cooling process. By determining even the distribution, it is possible to accurately predict the density distribution and size distribution of the oxygen precipitation nuclei in the single crystal. That is, by considering the effect of slow cooling and rapid cooling of the single crystal in the cooling process of the single crystal separated from the melt, the density distribution and size distribution of void defects consisting of voids and inner wall oxide films are obtained, and oxygen precipitation Calculate the amount of nuclei generated per unit time and this critical radius, and when the amount of oxygen precipitated nuclei generated per unit time exceeds zero, stop calculating the void radius and inner wall oxide thickness, The radius of the oxygen precipitation nuclei is obtained using a computer. As a result, the density distribution of oxygen precipitation nuclei in a single crystal can be accurately predicted by determining the amount of oxygen precipitation nuclei generated per unit time, and the size of the oxygen precipitation nuclei in the crystal can be determined by determining the radius of the oxygen precipitation nuclei. The distribution can be predicted accurately.
[Brief description of the drawings]
FIG. 1 is a flowchart showing a first stage of a simulation method of density distribution and size distribution of void defects in a silicon single crystal according to an embodiment of the present invention.
FIG. 2 is a flowchart showing a second stage of a method for simulating the density distribution and size distribution of the void defect.
FIG. 3 is a flowchart showing a third stage of the simulation method of the density distribution and size distribution of the void defect.
FIG. 4 is a flowchart showing a fourth stage of a simulation method of the density distribution and size distribution of the void defect.
FIG. 5 is a cross-sectional view of a main part of a silicon single crystal pulling machine having a mesh structure of the silicon melt according to the present invention.
FIG. 6 is a schematic view of a void and an inner wall oxide film in the silicon single crystal.
FIG. 7 is a cross-sectional view of a principal part of a silicon single crystal pulling machine having a mesh structure of a silicon melt of a conventional example.
FIG. 8 is a diagram showing the size distribution and density distribution of oxygen precipitation nuclei obtained using a computer when the radial position of the silicon single crystal is changed and when the pulling speed of the silicon single crystal is changed.
[Explanation of symbols]
11 Silicon single crystal pulling machine
12 Silicon melt
14 Silicon single crystal
21 void
22 Inner wall oxide film
S Triple point of silicon

Claims

The hot zone of the puller (11) is modeled with a mesh structure from the start of pulling of the single crystal (14) from the melt (12) by the puller (11) until the cooling of the single crystal (14) is completed. A first step to:
The mesh is grouped for each member of the hot zone and the pulling length of the single crystal (14) and the pulling speed of the single crystal (14) corresponding to the pulling length together with the physical property values of the members with respect to the gathered mesh A second step of entering each into the computer;
A third step of determining 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 on the assumption 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 fourth step of further obtaining an internal temperature distribution of the melt (12) considering convection by concatenating and solving a flow model equation and Navier-Stokes equations;
A fifth 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 third step to the fifth 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 sixth 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 in stages, and the temperature distribution in the pulling machine (11) is calculated repeatedly from the first step to the sixth step to calculate the single crystal (14 ) To obtain the coordinates and temperature of the mesh and input these data to the computer, respectively,
A 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) is divided at predetermined intervals, and a seventh step is performed for each divided time interval. An eighth step of determining the pulling length and pulling 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
Obtain the density distribution of the vacancies and interstitial atoms after the predetermined time interval by solving the diffusion equation based on the diffusion coefficients and boundary conditions of the vacancies and interstitial atoms in the single crystal (14) 9th step;
A tenth step of determining a void formation start temperature based on the density distribution of the pores;
An eleventh step of determining the density of the void (21) 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 (21);
A twelfth step of determining a radius of the void (21) 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 (21);
A thirteenth step of determining a formation start temperature of the inner wall oxide film (22) grown around the void (21);
The radius of the void (21) and the inner wall oxide film when the temperature at the lattice point of each mesh in the single crystal (14) gradually decreases and becomes lower than the formation start temperature of the inner wall oxide film (22). A fourteenth step for determining the thickness of (22) in relation to each other;
A fifteenth step of determining the amount of oxygen precipitation nuclei per unit time and the critical radius based on the density distribution of oxygen in the single crystal (14), the density distribution of the vacancies, and the temperature distribution;
When the generation amount per unit time of the oxygen precipitation nuclei exceeds zero, the calculation for determining the radius of the void (21) and the thickness of the inner wall oxide film (22) is stopped, and the radius of the oxygen precipitation nuclei is set. A 16th step to be obtained;
A method of simulating the density distribution and size distribution of oxygen precipitation nuclei in a single crystal using a computer including a seventeenth step in which steps 8 to 16 are repeated until the cooling of the single crystal (14) is completed.