JP4403722B2

JP4403722B2 - Simulation method for density distribution and size distribution of void defect in silicon single crystal

Info

Publication number: JP4403722B2
Application number: JP2003150262A
Authority: JP
Inventors: 浩之介北村; 純古川; 直樹小野
Original assignee: Sumco Corp
Current assignee: Sumco Corp
Priority date: 2003-05-28
Filing date: 2003-05-28
Publication date: 2010-01-27
Anticipated expiration: 2023-05-28
Also published as: JP2004356253A

Description

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

【００１４】
上記式（１）において、ρは各部材の密度、ｃは各部材の比熱、Ｔは各部材の各メッシュ点での絶対温度、ｔは時間である。またλ_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はボイド２１の形成開始温度である（図６）。更に式（９）及び（１０）のそれぞれ右側の第１項はフィックの拡散式であり、右側の第１項中のＤ_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はボイド２１の形成開始温度である。更に式（１３）の右側の第１項はフィックの拡散式であり、右側の第１項中のＤ_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と同一となる。
【００３９】
次に第１３ステップとして、ボイド２１の周囲に成長する内壁酸化膜２２の形成開始温度Ｔ_pを求める。内壁酸化膜２２の形成開始温度Ｔ_pがボイド２１の形成開始温度Ｔ_vより小さい場合には、酸素と空孔が結合し、ボイド２１表面に内壁酸化膜２２が成長する（図６）。本発明のモデルでは、上記内壁酸化膜２２はボイド２１が発生するとすぐに成長するものとして扱った。ここで、「すぐに」とは、「ボイドが発生した次の繰返し計算から」という意味である。従って、内壁酸化膜２２の形成開始温度Ｔ_pはボイド２１の形成開始温度Ｔ_vと同一となる。
【００４０】
上記式（１７）において、Ｄ_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＝ｘ／（ΩＳｉＯｘ）と定義する。このρ_pの定義式において、Ω_SiOxはＳｉＯｘの分子容量（molecular volume）であり、その単位は分子量／密度である。またｘはＳｉＯｘのｘに対応し、Ωは定数である。
【数１６】

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

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

【００４８】
上記式（２３）及び式（２４）において、Ω_SiOxはＳｉＯｘの分子容量（molecular volume）である。また式（２３）及び式（２４）において、内壁酸化膜の外周面での空孔平衡密度はＲ_pと無関係であり、Ｃ_veif＝Ｃ_veかつＣ_Oeif＝Ｃ_Oeであると想定している。式（２３）又は式（２４）から内壁酸化膜２２の外半径Ｒ_pを求め、内壁酸化膜２２の外半径Ｒ_pからボイド２１の半径ｒ_vを引くことにより、内壁酸化膜２２の厚さｄを求める。上記第８ステップから第１４ステップをシリコン単結晶１４の冷却が完了するまで繰返す（第１５ステップ）。なお、式（９）〜式（２４）は互いに関連させてコンピュータにより解く。
【００４９】
上述のように、シリコン融液１２の対流を考慮してシリコン融液１２から成長するシリコン単結晶１４内の温度分布を求めた後に、シリコン融液１２から切離されたシリコン単結晶１４の冷却過程、即ちシリコン単結晶１４がシリコン融液１２から切離された後のシリコン単結晶１４の引上げ速度を考慮し、シリコン単結晶１４の徐冷及び急冷の効果を結果に反映して解析することにより、シリコン単結晶１４内のボイド２１及び内壁酸化膜２２からなるボイド欠陥の密度分布及びサイズ分布を正確に予測できる。この結果、シリコン融液１２から引上げられるシリコン単結晶１４内の欠陥の密度分布及びサイズ分布を所望の分布にするために、引上げ機１１の設計段階で引上げ機１１の構造を検討することができる。
【００５０】
【実施例】
次に本発明の実施例を比較例とともに詳しく説明する。
＜実施例１＞
図５及び図６に示すように、石英るつぼ１３に貯留されたシリコン融液１２から直径６インチのシリコン単結晶１４を引上げる場合の、シリコン単結晶１４内のボイド２１及び内壁酸化膜２２からなるボイド欠陥の密度分布及びサイズ分布を、図１〜図４のフローチャートに基づくシミュレーション方法により求めた。
【００５１】
即ち、シリコン単結晶引上げ機１１のホットゾーンをメッシュ構造でモデル化した。ここで、シリコン融液１２のシリコン単結晶１４直下のシリコン単結晶１４の径方向のメッシュを０．７５ｍｍに設定し、シリコン融液１２のシリコン単結晶１４直下以外のシリコン単結晶１４の径方向のメッシュを１〜５ｍｍに設定した。またシリコン融液１２のシリコン単結晶１４の長手方向のメッシュを０．２５〜５ｍｍに設定し、乱流モデル式の乱流パラメータＣとして０．４５を用いた。更にシリコン単結晶１４の引上げ開始時ｔ₀から冷却完了時ｔ₁までの引上げ長及び引上げ高さの段階的な変更を５０ｍｍずつとした。このような条件下で、シリコン融液１２の対流を考慮してシリコン単結晶１４内の温度分布を求めた。
【００５２】
次にシリコン単結晶１４の引上げ開始時ｔ₀から冷却完了時ｔ₁までの時間を所定の間隔Δｔ秒（微小な時間の間隔）毎に区切り、この時間間隔Δｔ秒毎に上記シリコン単結晶１４のメッシュの座標及び温度のデータから、シリコン単結晶１４の引上げ長及び引上げ高さとシリコン単結晶内１４の温度分布を求め、更にボイド２１の形成開始温度を求めてボイド２１の密度分布及びサイズ分布を求めた。即ち、シリコン単結晶１４をシリコン融液１２から切離した後のシリコン単結晶１４の徐冷及び急冷を考慮して、シリコン単結晶１４内のボイド２１の密度分布及びサイズ分布をコンピュータを用いてそれぞれ求めた。
【００５３】
一方、ボイド２１の周囲に成長する内壁酸化膜２２は、ボイド２１が発生するとすぐに成長するものとして扱ったので、内壁酸化膜２２の形成開始温度Ｔ_pをボイド２１の形成開始温度Ｔ_vと同一とした。ここで、「すぐに」とは、「ボイドが発生した次の繰返し計算から」という意味である。次に内壁酸化膜２２の形成開始温度より低くなったときの内壁酸化膜２２の外半径Ｒ_pをコンピュータを用いて求めた。更に内壁酸化膜２２の外半径Ｒ_pをボイド２１の半径から引くことにより内壁酸化膜２２の厚さｄを求めた。なお、内壁酸化膜の密度分布は上記ボイドの密度分布と同一とした。上記ボイド及び内壁酸化膜からなるボイド欠陥のサイズ分布、即ち内壁酸化膜２２の外半径Ｒ_pのシリコン単結晶の半径方向への分布を、シリコン単結晶の引上げ速度を変えて算出した結果を図８（ａ）及び（ｂ）に実線で示す。また上記ボイド及び内壁酸化膜からなるボイド欠陥の密度分布のシリコン単結晶の半径方向への分布を、シリコン単結晶の引上げ速度を変えて算出した結果を図９（ａ）及び（ｂ）に実線で示す。
【００５４】
＜比較例１＞
実施例１と同一形状のシリコン単結晶を同一の条件で引上げたときのＬＳＴＤ（Laser Scattering Tomograph Defects）を測定した。このＬＳＴＤのサイズ分布を、シリコン単結晶の引上げ速度を変えて測定した結果を図８（ａ）及び（ｂ）に黒い点で示す。またＬＳＴＤの密度分布を、シリコン単結晶の引上げ速度を変えて測定した結果を図９（ａ）及び（ｂ）に黒い点で示す。なお、上記ＬＳＴＤは、空孔型点欠陥の凝集体の一種であり、シリコン単結晶内に赤外線を照射したときにシリコンとは異なる屈折率を有する散乱光を発する源をいう。なお、上記ＬＳＴＤのサイズ分布及び密度分布は、上記実施例１のボイド欠陥のサイズ分布及び密度分布に対応させてもよいと考える。またＬＳＴＤは検出限界があるため、所定値より小さいサイズのものは検出できなかった。
【００５５】
＜評価＞
図８（ａ）及び（ｂ）から明らかなように、実施例１のボイド欠陥のサイズ分布と比較例１のＬＳＴＤのサイズ分布とは、サイズの検出限界以上で略一致した。
また図９（ａ）及び（ｂ）から明らかなように、実施例１のボイド欠陥の密度分布と比較例１のＬＳＴＤの密度分布とは、範囲▲１▼において略一致したけれども、範囲▲２▼では一致しなかった。範囲▲２▼において密度分布が一致しなかったのは、範囲▲２▼では検出可能なボイドと検出不能なボイドが入り交じっているため、比較例１のＬＳＴＤの密度分布が実施例１のボイド欠陥の密度分布より少なくなってしまったものと考えられる。
【００５６】
【発明の効果】
以上述べたように、本発明によれば、シリコン融液の対流を考慮してシリコン融液から成長するシリコン単結晶内の温度分布をコンピュータを用いて求めるだけでなく、更に冷却過程におけるシリコン単結晶内の温度分布までも求めることによって、即ちシリコン融液から切離されたシリコン単結晶の冷却過程におけるシリコン単結晶の徐冷及び急冷の効果を考慮することによって、ボイドの密度分布をコンピュータを用いて求めた後に、ボイドの半径と、ボイドの周囲に成長する内壁酸化膜の厚さとを互いに関連させてコンピュータを用いて求める。この結果、シリコン単結晶内のボイド及び内壁酸化膜からなるボイド欠陥の密度分布及びサイズ分布を正確に予測できる。
【図面の簡単な説明】
【図１】本発明実施形態シリコン単結晶内のボイド欠陥の密度分布及びサイズ分布のシミュレーション方法の第１段を示すフローチャート。
【図２】そのボイド欠陥の密度分布及びサイズ分布のシミュレーション方法の第２段を示すフローチャート。
【図３】そのボイド欠陥の密度分布及びサイズ分布のシミュレーション方法の第３段を示すフローチャート。
【図４】そのボイド欠陥の密度分布及びサイズ分布のシミュレーション方法の第４段を示すフローチャート。
【図５】本発明のシリコン融液をメッシュ構造としたシリコン単結晶の引上げ機の要部断面図。
【図６】そのシリコン単結晶内のボイド及び内壁酸化膜の模式図。
【図７】従来例のシリコン融液をメッシュ構造としたシリコン単結晶の引上げ機の要部断面図。
【図８】実施例１のボイド欠陥のシリコン単結晶の半径方向へのサイズ分布を、引上げ速度を変えてコンピュータにより求めた値と、比較例１のＬＳＴＤのシリコン単結晶の半径方向へのサイズ分布を、引上げ速度を変えて実際に測定した値とを示す図。
【図９】実施例１のボイド欠陥のシリコン単結晶の半径方向への密度分布を、引上げ速度を変えてコンピュータにより求めた値と、比較例１のＬＳＴＤのシリコン単結晶の半径方向への密度分布を、引上げ速度を変えて実際に測定した値とを示す図。
【符号の説明】
１１シリコン単結晶引上げ機
１２シリコン融液
１４シリコン単結晶
２１ボイド
２２内壁酸化膜
Ｓシリコンの三重点[0001]
BACKGROUND OF THE INVENTION
  The present invention is a silico pulled up by the Czochralski (hereinafter referred to as CZ) method.SimpleThe present invention relates to a method for computer simulation of density distribution and size distribution of glow-in defects in a crystal. The glow-in defect referred to here is a void defect composed of a void and an inner wall oxide film. For example, a crystal-origin particle (hereinafter referred to as COP) observed on the wafer surface or an inside of the crystal. This corresponds to LSTD (Laser Scattering Tomograph Defects).
[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 coordinates and temperature of the mesh of the silicon single crystal 4 are obtained from the internal temperature distribution. A method is known in which the density distribution of the interstitial silicon and vacancies is determined by using a computer by solving a diffusion equation based on the diffusion coefficients 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, there is a problem that the density distribution of interstitial silicon and vacancies in the silicon single crystal is significantly different from the actually measured value. Further, in the above conventional simulation method, defects in the silicon single crystal generated at a relatively high temperature during the crystal cooling process, for example, types of voids and oxygen precipitates, density distribution and size distribution have not been 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 purpose of the present invention is tosiliconIn consideration of melt convectionsiliconAfter analyzing the temperature distribution in the single crystal,siliconDisconnected from the meltsiliconBy analyzing the cooling process of a single crystal,siliconIt is possible to accurately predict the density distribution and size distribution of void defects consisting of voids in the single crystal and inner wall oxide film.siliconThe object is to provide a method for simulating the density distribution and size distribution of void defects in a single crystal.
[0007]
[Means for Solving the Problems]
  The invention according to claim 1 is based on a pulling machine 11 as shown in FIGS.siliconSingle crystal 14siliconFrom the start of pulling up from the melt 12siliconThe first step of modeling the hot zone of the puller 11 until the completion of cooling of the single crystal 14 with a mesh structure, the mesh for each member of the hot zone, and the physical property values of each member for the combined meshsiliconCorresponding to the pulling length of the single crystal 14 and this pulling lengthsiliconPulling speed of single crystal 14And turbulent flow parameters of silicon melt 12A second step of inputting the temperature 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 emissivity of each member, and the surface temperature distribution and the thermal conductivity of each member. After finding the internal temperature distribution of each member by solving the heat conduction equationsiliconConvection was taken into account by concatenating and solving the turbulent model and Navier-Stokes equations obtained on the assumption that the melt 12 is turbulent.siliconA fourth step for further determining the internal temperature distribution of the melt 12;

silicon

Single crystal 14 andsiliconThe solid-liquid interface shape of the melt 12siliconThe fifth step, which is determined according to the isotherm including the triple point S of the single crystal, and the third step to the fifth step are the triple points S.siliconCalculate the temperature distribution in the puller 11 repeatedly until the melting point of the single crystal 14 is reached.siliconA sixth step of obtaining coordinates and temperature of the mesh of the single crystal 14 and inputting these data to a computer, respectively;siliconBy changing the pulling length and pulling height of the single crystal 14 step by step, the temperature distribution in the pulling machine 11 is calculated repeatedly from the first step to the sixth step.siliconA seventh step of obtaining coordinates and temperature of the mesh of the single crystal 14 and inputting these data to a computer, respectively;siliconSingle crystal 14siliconFrom the start of pulling up from the melt 12siliconThe time until the completion of the cooling of the single crystal 14 was divided at predetermined intervals, and found at the seventh step for each divided time interval.siliconFrom the coordinates and temperature data of the mesh of the single crystal 14siliconThe pulling length and pulling height of the single crystal 14siliconAn eighth step for obtaining a temperature distribution in the single crystal 14;siliconA ninth step for determining a density distribution of vacancies and interstitial atoms after a predetermined time interval by solving a diffusion equation based on diffusion coefficients and boundary conditions of vacancies and interstitial atoms in the single crystal 14; A tenth step of determining the formation start temperature of the void 21 based on the density distribution of the holes;siliconAn eleventh step of obtaining the density of the voids 21 when the temperature at the lattice points of the respective meshes in the single crystal 14 gradually decreases to the formation start temperature of the voids 21;siliconA twelfth step of obtaining the radius of the void 21 when the temperature at the lattice point of each mesh in the single crystal 14 is lower than the formation start temperature of the void 21, and the formation start temperature of the inner wall oxide film 22 grown around the void 21 13th step for obtainingsiliconThe radius of the void 21 and the thickness of the inner wall oxide film 22 are obtained in relation to each other 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. 14th step and 8th to 14th stepssiliconUsing a computer, the fifteenth step being repeated until cooling of the single crystal 14 is completedsiliconThis is a method for simulating the density distribution and size distribution of void defects in a single crystal.
[0008]
In the simulation method according to claim 1,siliconConsidering the convection of the melt 12siliconGrows from melt 12siliconIn addition to obtaining the temperature distribution in the single crystal 14, further in the cooling processsiliconBy determining even the temperature distribution in the single crystal 14,siliconSeparated from the melt 12siliconIn the cooling process of the single crystal 14siliconBy analyzing the effect of slow cooling and rapid cooling of the single crystal 14,siliconIt is possible to accurately predict the density distribution and size distribution of the void defects formed of the void 21 and the inner wall oxide film 22 in the single crystal 14.
[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 for simulating the density distribution and size distribution of void defects 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 pulled up to (for example, 100 mm), that is, a chamber, a quartz crucible 13, a silicon melt 12, a silicon single crystal 14, a graphite susceptor, a 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 wall of the quartz crucible 13 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 in the inside is calculated to determine the coordinates and temperature of the mesh of the silicon single crystal, 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 level 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 Δt seconds has been 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). Further, 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 formula (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 consumption rate of vacancies per oxygen atom, and D_OIs the diffusion coefficient of oxygen, C_OIs the oxygen density.
[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. Further, 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 is D_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 diameter of the void 21 and this critical diameter 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]
Next, as a thirteenth step, a formation start temperature T of the inner wall oxide film 22 growing around the void 21 is obtained._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 is generated.” Therefore, the formation start temperature T of the inner wall oxide film 22 is determined._pIs the formation start temperature T of the void 21_vIs the same as
[0039]
Next, as a thirteenth step, a formation start temperature T of the inner wall oxide film 22 growing around the void 21 is obtained._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 that a void has occurred."It means that. Accordingly, the formation start temperature T of the inner wall oxide film 22 is increased._pIs the formation start temperature T of the void 21_vIs the same as
[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.
[0041]
In addition, the voids are formed from outside the inner wall oxide film 22 to the outer peripheral surface (radius R_pFlux flowing into the spherical surface_vIs obtained from the following equation (18).
[Formula 13]

[0042]
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]

[0043]
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.
[0044]
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 distortion-free inner wall oxide film 22 cannot be grown, so that the voids flow into the oxide film from inside the void 21 and the outer diameter of the void 21 becomes small. 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]

[0045]
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 the oxide film._OTherefore, the following equation (21) is established. In equation (21), ρ_p= X / (ΩSiOx). This ρ_pIn the definition formula, Ω_SiOxIs the molecular volume of SiOx, the unit of which is molecular weight / density. Further, x corresponds to x of SiOx, and Ω is a constant.
[Expression 16]

[0046]
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]

[0047]
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 amount of vacancies or oxygen consumed 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]

[0048]
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. The eighth step to the fourteenth step are repeated until the cooling of the silicon single crystal 14 is completed (fifteenth step). Equations (9) to (24) are solved by a computer in association with each other.
[0049]
  As described above, after obtaining 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. In consideration of the process, that is, the pulling speed of the silicon single crystal 14 after the silicon single crystal 14 is separated from the silicon melt 12, the effect of slow cooling and rapid cooling of the silicon single crystal 14 is reflected in the analysis. Thus, it is possible to accurately predict the density distribution and size distribution of the void defects formed of the voids 21 and the inner wall oxide film 22 in the silicon single crystal 14. As a result, the structure of the puller 11 can be studied at the design stage of the puller 11 in order to obtain the desired density distribution and size distribution of defects in the silicon single crystal 14 pulled from the silicon melt 12.The
[0050]
【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.
[0051]
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.
[0052]
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.
[0053]
  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 that a void has occurred."It means that. 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. Further, the outer radius R of the inner 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. Size distribution of void defects composed of the void and the inner wall oxide film, that is, the outer radius R of the inner wall oxide film 22._pThe results of calculating the distribution of the silicon single crystal in the radial direction by changing the pulling speed of the silicon single crystal are shown by solid lines in FIGS. Also, the results of calculating the distribution of the void defect density distribution composed of the voids and the inner wall oxide film in the radial direction of the silicon single crystal by changing the pulling speed of the silicon single crystal are shown by solid lines in FIGS. 9A and 9B. Indicated by
[0054]
<Comparative Example 1>
LSTD (Laser Scattering Tomograph Defects) when a silicon single crystal having the same shape as in Example 1 was pulled under the same conditions was measured. The result of measuring the size distribution of this LSTD while changing the pulling rate of the silicon single crystal is shown by black dots in FIGS. 8 (a) and 8 (b). Further, the results of measuring the density distribution of LSTD by changing the pulling rate of the silicon single crystal are shown by black dots in FIGS. 9 (a) and 9 (b). The LSTD is a kind of agglomerates of vacancy-type point defects and refers to a source that emits scattered light having a refractive index different from that of silicon when an infrared ray is irradiated into a silicon single crystal. Note that the size distribution and density distribution of the LSTD may correspond to the size distribution and density distribution of the void defects of the first embodiment. In addition, since LSTD has a detection limit, a size smaller than a predetermined value could not be detected.
[0055]
<Evaluation>
As is apparent from FIGS. 8A and 8B, the size distribution of the void defect in Example 1 and the size distribution of the LSTD in Comparative Example 1 substantially coincided with each other above the size detection limit.
9A and 9B, the void defect density distribution of Example 1 and the LSTD density distribution of Comparative Example 1 substantially coincide with each other in the range {circle around (1)}, but the range {circle around (2)}. ▼ did not agree. In the range (2), the density distributions did not coincide with each other because the detectable voids and the non-detectable voids were mixed in the range (2). Therefore, the density distribution of the LSTD in Comparative Example 1 was the void in Example 1. This is considered to be less than the defect density distribution.
[0056]
【The invention's effect】
  As described above, according to the present invention,siliconConsidering the convection of the meltsiliconGrowing from the meltsiliconNot only can the temperature distribution in a single crystal be determined using a computer,siliconBy determining even the temperature distribution in the single crystal,siliconDisconnected from the meltsiliconIn the cooling process of single crystalssiliconBy taking into account the effects of slow cooling and rapid cooling of the single crystal, the density distribution of the void is obtained using a computer, and then the void radius and the thickness of the inner wall oxide film grown around the void are correlated with each other. Find using a computer. As a result,siliconIt is possible to accurately predict the density distribution and size distribution of void defects formed of voids in the single crystal and the inner wall oxide film.
[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 a 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 diagram of voids 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 according to a conventional example.
FIG. 8 shows the size distribution in the radial direction of the silicon single crystal of the void defect of Example 1 obtained by a computer while changing the pulling rate, and the size of the LSTD of Comparative Example 1 in the radial direction of the silicon single crystal. The figure which shows distribution and the value actually measured by changing pulling speed.
9 shows the density distribution in the radial direction of the silicon single crystal of the void defect in Example 1 obtained by a computer while changing the pulling rate, and the density in the radial direction of the silicon single crystal of the LSTD in Comparative Example 1. FIG. The figure which shows distribution and the value actually measured by changing pulling speed.
[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) from the start of pulling of the silicon single crystal (14) from the silicon melt (12) by the puller (11) to the completion of cooling of the silicon single crystal (14) has a mesh structure The first step of modeling with
The summarized mesh for each member of the hot zone and the pulling length and the silicon single crystal corresponding to the pull-up length of the silicon single crystal with physical properties of each member for the gathered mesh (14) (14) A second step of inputting the pulling speed and the turbulent flow parameter of the silicon melt (12) to a 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;
Obtained by assuming that the silicon 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 silicon melt (12) considering convection by connecting and solving a turbulent model equation and Navier-Stokes equations;
A fifth step of obtaining a solid-liquid interface shape of the silicon single crystal (14) and the silicon melt (12) according to an isotherm including a triple point (S) of the silicon single crystal;
The third step to the fifth step are repeated until the triple point (S) reaches the melting point of the silicon single crystal (14), and the temperature distribution in the puller (11) is calculated to calculate the silicon single crystal (14 6) to obtain the coordinates and temperature of the mesh) and input these data to the computer,
The silicon single crystal (14) is gradually changed in pulling height and pulling height, the first step to the sixth step are repeated, and the temperature distribution in the pulling machine (11) is calculated to calculate the silicon single crystal. (14) obtaining the coordinates and temperature of the mesh, and inputting these data to the computer, respectively,
The time from the start of pulling of the silicon single crystal (14) from the silicon melt (12) to the completion of cooling of the silicon single crystal (14) is divided at predetermined intervals, and the time intervals are divided. obtaining a temperature distribution from the mesh coordinates and temperature data of the pulling length and the pulling height and the silicon single crystal (14) of the silicon single crystal (14) the silicon single crystal (14) obtained in the seventh step The eighth step;
The density distribution of vacancies and interstitial atoms after elapse of the predetermined time interval is solved by solving a diffusion equation based on diffusion coefficients and boundary conditions of vacancies and interstitial atoms in the silicon single crystal (14). The ninth step to be determined;
A tenth step of determining the formation start temperature of the void (21) based on the density distribution of the holes;
An eleventh step of determining the density of the void (21) when the temperature at the lattice point of each mesh in the silicon single crystal (14) gradually decreases to the formation start temperature of the void (21);
A twelfth step of obtaining a radius of the void (21) when a temperature at a lattice point of each mesh in the silicon 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 oxidation when the temperature at the lattice point of each mesh in the silicon single crystal (14) is gradually decreased to be lower than the formation start temperature of the inner wall oxide film (22). A fourteenth step for determining the thickness of the film (22) in relation to each other;
A method of simulating the density distribution and size distribution of void defects in a silicon single crystal using a computer comprising: a fifteenth step in which steps 8 to 14 are repeated until the cooling of the silicon single crystal (14) is completed.