JP4192704B2 - A simulation method to maximize the defect-free region of a single crystal - Google Patents

Description

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

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

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

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

ここで、Ｋ₁及びＫ₂は定数であり、Ｅ_i及びＥ_vはそれぞれ格子間シリコン及び空孔の形成エネルギーであり、Ｃ_v ^e及びＣ_i ^eの肩付き文字ｅは平衡量、ｋはボルツマン定数、Ｔは絶対温度を意味する。
上記平衡式は時間で微分され、空孔及び格子間シリコンに対してそれぞれ次の式（９）及び式（１０）になる。
【００３１】
【数６】

ここで、Θ(ｘ)はヘビサイド関数（Heaviside function）である。即ち、ｘ＜０のときΘ(ｘ)＝０であり、かつｘ＞０のときΘ(ｘ)＝１である。またＴ_nは高温酸素析出物の形成開始温度Ｔ_pとボイドの形成開始温度Ｔ_v0とを比較したときの高い方の温度である。更に式（９）及び（１０）のそれぞれ右側の第１項はフィックの拡散式であり、右側の第１項中のＤ_v及びＤ_iは、次の式（１１）及び（１２）に表される拡散係数である。
【００３２】
【数７】

ここで△Ｅ_v及び△Ｅ_iはそれぞれ空孔及び格子間シリコンの活性化エネルギーであり、ｄ_v及びｄ_iはそれぞれ定数である。また式（９）及び式（１０）のそれぞれ右側の第２項中の
【数８】

は熱拡散による空孔及び格子間シリコンの活性化エネルギーであり、ｋ_Bはボルツマン定数である。式（９）及び式（１０）のそれぞれ右側の第３項のｋ_ivは空孔及び格子間シリコンペアの再結合定数である。式（９）及び式（１０）のそれぞれ右側の第４項のＮ_v0はボイドの濃度であり、ｒ_v0はボイドの半径であり、更に式（９）の右側の第５項のＮ_pは高温酸素析出物の濃度であり、Ｒ_pは高温酸素析出物の半径であり、γはＳｉＯ₂析出物が歪まずに析出するために必要な酸素１原子当りの空孔消費量である。
【００３３】
一方、上記式（９）が成立つのは、空孔が析出するための流速が十分大きく、ＳｉマトリックスとＳｉＯ₂との単位質量当りの体積差を埋められる場合、即ちＤ_v（Ｃ_v−Ｃ_v ^e）≧γＤ₀Ｃ₀の場合である。上記以外の場合には、次式（１３）が成立つ。
【００３４】
【数９】

次に第１１ステップとして、上記拡散方程式を解くことにより求めた空孔の濃度Ｃ_v分布に基づいて、ボイドの形成開始温度Ｔ_v0を次の式（１４）から求める。
【００３５】
【数１０】

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

更に第１４ステップとして、シリコン単結晶１４内のそれぞれのメッシュの格子点における温度がボイドの形成開始温度Ｔ_v0より低いときのボイドの半径ｒ_v0を、次の式（１６）から求める。
【００３８】
【数１２】

ここで、ｔ₁はシリコン単結晶１４のメッシュの格子点における温度がボイドの形成開始温度Ｔ_v0まで低下したときの時刻であり、ｒ_crはボイドの臨界径である。上記第９〜第１１ステップ、第１３及び第１４ステップをシリコン単結晶１４が８００〜１０００℃の間の特定値、例えば９００℃以下に冷却するまで繰返す（第１７ステップ）。なお、式（９）〜式（１６）は連成してコンピュータにより解く。また第１５及び第１６ステップは後述する。
【００３９】
▲２▼ シリコン単結晶１４内の高温酸素析出物の濃度分布及びサイズ分布
第１２ステップに戻って、上記拡散方程式（９）及び（１０）を解くことにより求めた空孔の濃度Ｃ_v分布に基づいて、高温酸素析出物の形成開始温度Ｔ_pを次の式（１７）から計算により求める。
【００４０】
【数１３】

ここで、Ｃ₀は酸素濃度であり、Ｃ_0mはシリコン融液１２の融点Ｔ_mでの酸素平衡濃度であり、Ｅ₀は酸素溶解エネルギーである。またＥ_vは空孔形成エネルギーであり、σ_pはシリコン単結晶１４内のＳｉとＳｉＯ₂との界面エネルギーである。更にγはＳｉＯ₂析出物が歪まずに析出するために必要な酸素１原子当りの空孔消費量であり、その値は０．６８である。
【００４１】
次に第１５ステップとして、シリコン単結晶１４内のそれぞれのメッシュの格子点における温度が次第に低下して高温酸素析出物の形成開始温度Ｔ_pになったときに、次の近似式（１８）を用いて高温酸素析出物の濃度Ｎ_pを求める。なお、式（１８）において、ａ₂は定数である。
【００４２】
【数１４】

更に第１６ステップとして、シリコン単結晶１４内のそれぞれのメッシュの格子点における温度が高温酸素析出物の形成開始温度Ｔ_pより低いときの高温酸素析出物の半径Ｒ_pを、次の式（１９）から求める。
【００４３】
【数１５】

ここで、ｔ₂はシリコン単結晶１４のメッシュの格子点における温度が高温酸素析出物の形成開始温度Ｔ_pまで低下したときの時刻であり、Ｒ_crは高温酸素析出物の臨界径である。
【００４４】
一方、上記式（１９）が成立つのは、空孔が析出するための流速が十分大きく、ＳｉマトリックスとＳｉＯ₂との単位質量当りの体積差を埋められる場合、即ちＤ_v（Ｃ_v−Ｃ_v ^e）≧γＤ₀Ｃ₀の場合である。上記以外の場合には、次式（２０）が成立つ。
【００４５】
【数１６】

上記第９、第１０、第１２、第１５及び第１６ステップをシリコン単結晶１４が８００〜１０００℃の間の特定値、例えば９００℃以下に冷却するまで繰返す（第１７ステップ）。なお、上記式（９）〜式（１３）及び式（１７）〜式（２０）は連成してコンピュータにより解く。
【００４６】
第１８ステップとして、上記高温酸素析出物の存在するＰバンドとボイド及び高温酸素析出物のいずれも存在しない無欠陥領域とを区画し、かつ無欠陥領域がＰバンドより単結晶ボトム側に位置するように区画する第１等濃度線ＨＣ１_Xを計算により求める（図７）。次に第１９ステップとして、高濃度酸素析出物が存在するＢバンドと無欠陥領域とを区画し、かつ無欠陥領域がＢバンドより単結晶トップ側に位置するように区画する第１分布線ＢＣ１_Xを計算により求める。上記第１分布線は、格子間シリコン濃度がシリコン融点での格子間シリコンの平衡濃度に対して０．１２〜０．１３％の範囲の特定値、例えば０．１２６％である等濃度線をいう。ここで格子間シリコン濃度０．１２〜０．１３％の範囲の特定値、例えば０．１２６％とは、シミュレーションを行って得られた格子間シリコン濃度のうちライフタイムのマップ等から観察されるＢバンドに対応する濃度である。またＢバンドとは、格子間シリコンの凝集体が核となり熱処理によって酸素析出が高濃度に発生している領域をいう。
【００４７】
その後、第２０ステップとして、第１等濃度線ＨＣ１_Xの変曲点の最大値と第１分布線ＢＣ１_Xの変曲点の最小値との差ΔＺ₁(X)を求める。第１等濃度線ＨＣ１_Xの変曲点が例えば３つある場合には、各変曲点Ｑ１(1,X)、Ｑ１(2,X)及びＱ１(3,X)の座標を(ｒ_Q1(1,X)，ｚ_Q1(1,X))、(ｒ_Q1(2,X)，ｚ_Q1(2,X))及び(ｒ_Q1(3,X)，ｚ_Q1(3,X))とすると、これらの座標のうちｚ座標ｚ_Q1(1,X)、ｚ_Q1(2,X)及びｚ_Q1(3,X)のうちの最大値ｚ_Q1(1,X)を求める。また第１分布線ＢＣ１_Xの変曲点が例えば５つある場合には、各変曲点Ｓ１(1,X)、Ｓ１(2,X)、Ｓ１(3,X)、Ｓ１(4,X)及びＳ１(5,X)の座標を(ｒ_S1(1,X)，ｚ_S1(1,X))、(ｒ_S1(2,X)，ｚ_S1(2,X))、(ｒ_S1(3,X)，ｚ_S1(3,X))、(ｒ_S1(4,X)，ｚ_S1(4,X))及び(ｒ_S1(5,X)，ｚ_S1(5,X))とすると、これらの座標のうちｚ座標ｚ_S1(1,X)、ｚ_S1(2,X)、ｚ_S1(3,X)、ｚ_S1(4,X)及びｚ_S1(5,X)のうちの最小値ｚ_S1(3,X)を求める。そして第１等濃度線ＨＣ１_Xの変曲点の最大値ｚ_Q1(1,X)と第１分布線ＢＣ１_Xの変曲点の最小値ｚ_S1(3,X)との差ΔＺ₁(X)を求める。
【００４８】
第２１ステップとして、単結晶製造条件のパラメータをＰ₂に変え第２ステップから第２０ステップを実行し、単結晶製造条件のパラメータをＰ₃に変え第２ステップから第２０ステップを実行するというように、単結晶製造条件のパラメータがＰ_Nになるまで上記第２ステップから第２０ステップを繰返した後に、第１等濃度線ＨＣ１_Xの変曲点の最大値と第１分布線ＢＣ１_Xの変曲点の最小値との差ΔＺ₁(X)を求め、更にこの差ΔＺ₁(X)が最も大きくなる単結晶製造条件を計算により求める。
【００４９】
このように計算して求められたシリコン単結晶引上げ機１１における単結晶製造条件は、シリコン単結晶１４の引上げ方向及び半径方向に無欠陥領域が最も拡大する条件であり、この単結晶製造条件を上記計算により正確に求めることができる。この結果、上述のようにコンピュータを用いた計算によりシリコン単結晶１４の無欠陥領域が最も拡大する単結晶製造条件を設定することができるので、実際にシリコン単結晶を引上げてライフタイム等を測定する実験は上記単結晶製造条件が無欠陥領域が実際に拡大しているか否かの確認のためにだけ行えば済む。従って、上記実験の回数を低減できるので、比較的少ない時間と労力で、最適な単結晶製造条件を得ることができる。なお、図７から明らかなように、上記差ΔＺ₁(X)が大きくなると、シリコン単結晶１４の引上げ方向及び半径方向にわたって無欠陥領域とするためのシリコン単結晶１４の引上げ速度の幅、即ちピュアマージンが大きくなるので、全長のわたって無欠陥のシリコン単結晶１４を比較的容易に引上げることができる。
【００５０】
次にシリコン単結晶引上げ機１１におけるシリコン単結晶１４の無欠陥領域を最大化する第２のシミュレーション方法を図６〜図１２に基づいて説明する。
第１〜第１７ステップは、上記第１のシミュレーション方法の第１〜第１７ステップと同一である。
【００５１】
第１８ステップとして、上記高温酸素析出物の存在するＰバンドとボイド及び高温酸素析出物のいずれも存在しない無欠陥領域とを区画し、かつ無欠陥領域がＰバンドより単結晶トップ側に位置するように区画する第２等濃度線ＨＣ２_Xを計算により求める（図７）。次に第１９ステップとして、高濃度酸素析出物が存在するＢバンドと無欠陥領域とを区画し、かつ無欠陥領域がＢバンドより単結晶ボトム側に位置するように区画する第２分布線ＢＣ２_Xを計算により求める。上記第２分布線は、格子間シリコン濃度がシリコン融点での格子間シリコンの平衡濃度に対して０．１２〜０．１３％の範囲の特定値、例えば０．１２６％である等濃度線をいう。ここで格子間シリコン濃度０．１２〜０．１３％の範囲の特定値、例えば０．１２６％とは、シミュレーションを行って得られた格子間シリコン濃度のうちライフタイムのマップ等から観察されるＢバンドに対応する濃度である。
【００５２】
その後、第２０ステップとして、第２等濃度線ＨＣ２_Xの変曲点の最小値と第２分布線ＢＣ２_Xの変曲点の最大値との差ΔＺ₂(X)を求める。第２等濃度線ＨＣ２_Xの変曲点が例えば３つある場合には、各変曲点Ｑ２(1,X)、Ｑ２(2,X)及びＱ２(3,X)の座標を(ｒ_Q2(1,X)，ｚ_Q2(1,X))、(ｒ_Q2(2,X)，ｚ_Q2(2,X))及び(ｒ_Q2(3,X)，ｚ_Q2(3,X))とすると、これらの座標のうちｚ座標ｚ_Q2(1,X)、ｚ_Q2(2,X)及びｚ_Q2(3,X)のうちの最小値ｚ_Q2(1,X)を求める。また第２分布線ＢＣ２_Xの変曲点が例えば５つある場合には、各変曲点Ｓ２(1,X)、Ｓ２(2,X)、Ｓ２(3,X)、Ｓ２(4,X)及びＳ２(5,X)の座標を(ｒ_S2(1,X)，ｚ_S2(1,X))、(ｒ_S2(2,X)，ｚ_S2(2,X))、(ｒ_S2(3,X)，ｚ_S2(3,X))、(ｒ_S2(4,X)，ｚ_S2(4,X))及び(ｒ_S2(5,X)，ｚ_S2(5,X))とすると、これらの座標のうちｚ座標ｚ_S2(1,X)、ｚ_S2(2,X)、ｚ_S2(3,X)、ｚ_S2(4,X)及びｚ_S2(5,X)のうちの最大値ｚ_S2(3,X)を求める。そして第２等濃度線ＨＣ２_Xの変曲点の最小値ｚ_Q2(1,X)と第２分布線ＢＣ２_Xの変曲点の最大値ｚ_S2(3,X)との差ΔＺ₂(X)を求める。
【００５３】
第２１ステップとして、単結晶製造条件のパラメータをＰ₂に変え第２ステップから第２０ステップを実行し、単結晶製造条件のパラメータをＰ₃に変え第２ステップから第２０ステップを実行するというように、単結晶製造条件のパラメータがＰ_Nになるまで上記第２ステップから第２０ステップを繰返した後に、第２等濃度線ＨＣ２_Xの変曲点の最小値と第２分布線ＢＣ２_Xの変曲点の最大値との差ΔＺ₂(X)を求め、更にこの差ΔＺ₂(X)が最も大きくなる単結晶製造条件を計算により求める。
【００５４】
このように計算して求められたシリコン単結晶引上げ機１１における単結晶製造条件は、シリコン単結晶１４の引上げ方向及び半径方向に無欠陥領域が最も拡大する条件であり、この単結晶製造条件を上記計算により正確に求めることができる。この結果、上述のようにコンピュータを用いた計算によりシリコン単結晶１４の無欠陥領域が最も拡大する単結晶製造条件を設定することができるので、実際にシリコン単結晶を引上げてライフタイム等を測定する実験は上記単結晶製造条件が無欠陥領域が実際に拡大しているか否かの確認のためにだけ行えば済む。従って、上記実験の回数を低減できるので、比較的少ない時間と労力で、最適な単結晶製造条件を得ることができる。なお、図７から明らかなように、上記差ΔＺ₂(X)が大きくなると、シリコン単結晶１４の引上げ方向及び半径方向にわたって無欠陥領域とするためのシリコン単結晶１４の引上げ速度の幅、即ちピュアマージンが大きくなるので、全長のわたって無欠陥のシリコン単結晶１４を比較的容易に引上げることができる。
なお、上記実施の形態では、単結晶としてシリコン単結晶を挙げたが、ＧａＡｓ単結晶，ＩｎＰ単結晶，ＺｎＳ単結晶若しくはＺｎＳｅ単結晶でもよい。
【００５５】
【発明の効果】
以上述べたように、本発明によれば、先ず所定の単結晶製造条件で融液の対流を考慮することにより、計算した結晶及び融液の固液界面形状が実際の単結晶の引上げ時の形状とほぼ一致するようにして単結晶の内部温度を求めるだけでなく、更に冷却過程における単結晶内の温度分布までも求めることによって、即ち融液から切離された単結晶の冷却過程における単結晶の徐冷及び急冷の効果を考慮することによって、単結晶内の欠陥の濃度分布及びサイズ分布をコンピュータを用いて求める。次に単結晶内の第１等濃度線と第１分布線を計算により求めて、第１等濃度線の変曲点の最大値と第１分布線の変曲点の最小値との差を計算により求めた後に、上記単結晶製造条件のパラメータを変え、上記と同様にして第１等濃度線の変曲点の最大値と第１分布線の変曲点の最小値との差を計算により求める。更に上記第１等濃度線の変曲点の最大値と第１分布線の変曲点の最小値との差が最も大きくなる単結晶製造条件を求める。これにより単結晶の引上げ方向及び半径方向に無欠陥領域が最も拡大する単結晶製造条件を計算により正確に求めることができる。このため上記無欠陥領域が最も拡大する単結晶製造条件を設定するのに数多くの実験を繰返す必要のあった従来の無欠陥の単結晶インゴットと比較して、本発明では、上記無欠陥領域が最も拡大する単結晶製造条件を設定するための実験の回数を低減できる。この結果、本発明では、比較的少ない時間と労力で、無欠陥領域が最も拡大する単結晶製造条件を得ることができる。
【００５６】
また単結晶内の第２等濃度線と第２分布線を計算により求めて、第２等濃度線の変曲点の最大値と第２分布線の変曲点の最小値との差を計算により求めた後に、上記単結晶製造条件のパラメータを変え、上記と同様にして第２等濃度線の変曲点の最大値と第２分布線の変曲点の最小値との差を計算により求め、更に上記第１等濃度線の変曲点の最大値と第１分布線の変曲点の最小値との差が最も大きくなる単結晶製造条件を求めても、上記と同様の効果が得られる。
【図面の簡単な説明】
【図１】本発明実施形態の第１のシミュレーション方法の第１段を示すフローチャート。
【図２】そのシミュレーション方法の第２段を示すフローチャート。
【図３】そのシミュレーション方法の第３段を示すフローチャート。
【図４】そのシミュレーション方法の第４段を示すフローチャート。
【図５】そのシミュレーション方法の第５段を示すフローチャート。
【図６】シリコン融液をメッシュ構造としたシリコン単結晶の引上げ機の要部断面図。
【図７】そのシリコン単結晶の引上げ速度を変化させたときのシリコン単結晶内の格子間シリコン及び空孔の分布を示す説明図。
【図８】本発明実施形態の第２のシミュレーション方法の第１段を示すフローチャート。
【図９】そのシミュレーション方法の第２段を示すフローチャート。
【図１０】そのシミュレーション方法の第３段を示すフローチャート。
【図１１】そのシミュレーション方法の第４段を示すフローチャート。
【図１２】そのシミュレーション方法の第５段を示すフローチャート。
【符号の説明】
１１ シリコン単結晶引上げ機
１２ シリコン融液
１４ シリコン単結晶
Ｓ シリコンの三重点[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a method for maximizing a defect-free region of a single crystal such as silicon pulled by a Czochralski (hereinafter referred to as CZ) method by computer simulation.
[0002]
[Prior art]
Factors that reduce device yields due to the recent miniaturization of semiconductor integrated circuits include crystal-origin particles (hereinafter referred to as COP), oxidation-induced stacking faults (hereinafter referred to as COP), The presence of microdefects in oxygen precipitates, which are the core of OISF, or interstitial-type large dislocation (hereinafter referred to as L / D).
[0003]
COP is a pit that appears on the wafer surface when a mirror-polished silicon wafer is washed with a mixed solution of ammonia and hydrogen peroxide. When this wafer is measured with a particle counter, the pits are detected as particles. This pit is caused by crystals. COP causes deterioration of electrical characteristics, for example, dielectric breakdown characteristics (Time Dependent dielectric Breakdown, TDDB) of oxide film, breakdown voltage characteristics of oxide film (Time Zero Dielectric Breakdown, TZDB), and the like. Further, if COP exists on the wafer surface, a step is generated in the device wiring process, which may cause disconnection. In addition, the element isolation portion also causes leakage and the like, thereby reducing the product yield.
[0004]
OISF is a defect that is manifested in a thermal oxidation process or the like when a semiconductor device is manufactured by introducing minute oxygen precipitation nuclei into a silicon single crystal during crystal growth. This OISF causes a failure such as increasing the leakage current of the device. L / D is also referred to as a dislocation cluster, or it is also referred to as a dislocation pit because a pit is generated when a silicon wafer having such a defect is immersed in a selective etching solution containing hydrofluoric acid as a main component. This L / D also causes deterioration of electrical characteristics such as leakage characteristics and isolation characteristics. From the above, it is necessary to reduce COP, OISF, and L / D from a silicon wafer used for manufacturing a semiconductor integrated circuit.
[0005]
A defect-free single crystal ingot having no COP, OISF, and L / D and a silicon wafer sliced from the single crystal ingot are disclosed (for example, see Patent Document 1). This defect-free single crystal ingot has a perfect region [P] when a perfect region in which no agglomerates of vacancy-type point defects and agglomerates of interstitial silicon type point defects exist in the single crystal ingot, respectively. P] is a single crystal ingot. The perfect region [P] includes a region [V] having point defects in which vacancy-type point defects are dominant and agglomerated in a single crystal ingot, and a point defect in which interstitial silicon type point defects are dominant and agglomerated. It is interposed between the region [I].
[0006]
[Patent Document 1]
Japanese Patent Laid-Open No. 11-1393 corresponding to US Pat. No. 6,045,610
[0007]
[Problems to be solved by the invention]
However, the conventional defect-free single crystal ingot disclosed in Japanese Patent Application Laid-Open No. 11-1393 has a void-type point defect aggregate and an interstitial silicon-type point defect aggregate over the entire length of the single crystal ingot. In order to set the pulling conditions that will not be performed, it is necessary to repeat a number of experiments in which the single crystal ingot is actually pulled and the lifetime measured for each type of puller and the thickness of the single crystal ingot. There was a problem that required a great deal of time and effort to set up.
In addition, in the conventional defect-free single crystal ingot, the V / G theory used for obtaining the point defect concentration, that is, the pulling rate of the single crystal is V (mm / min), and the solid crystals and the melt are solidified. When the temperature gradient near the liquid interface is G (° C./mm), V / G (mm²Because it is a one-dimensional model (model only in the pulling direction of the single crystal) that simplifies the Boronkov theory that controls (/ min · ° C), the point defect diffusion (outward diffusion) in the radial direction of the single crystal is considered. However, there is a problem that the concentration gradient of point defects in the radial direction of the single crystal cannot be obtained.
Further, in the conventional defect-free single crystal ingot, when the concentration of point defects is determined by V / G theory, it is necessary to model and calculate the convection of the melt in order to determine G with high accuracy. Nevertheless, since such consideration is not made at all, there is a problem that the concentration distribution of point defects in the single crystal cannot be obtained with high accuracy.
[0008]
The object of the present invention is to reduce the number of experiments of actually pulling up a single crystal, and thereby the defect-free region is expanded most in the pulling direction and radial direction of the single crystal with relatively little time and effort. An object of the present invention is to provide a simulation method capable of accurately obtaining single crystal production conditions by calculation.
[0009]
[Means for Solving the Problems]
As shown in FIGS. 1 to 7, the invention according to claim 1 is a parameter group P in which variables fed back as a single crystal production condition when pulling up the single crystal 14 by the puller 11 are varied at regular intervals.₁, P₂, ..., P_NAnd a parameter group P₁, P₂, ..., P_NParameter P₁A second step of modeling the hot zone of the puller 11 with a mesh structure from the start of pulling of the single crystal 14 from the melt to the completion of cooling of the single crystal 14 under the single crystal manufacturing conditions, and each member of the hot zone The parameter group P is set together with the physical property values of each member for the meshes that are grouped together.₁, P₂, ..., P_NParameter P of₁The third step of inputting the temperature into the computer, the fourth 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 heat based on the surface temperature distribution and thermal conductivity of each member. After obtaining the internal temperature distribution of each member by solving the conduction equation, the turbulence model equation obtained by assuming that the melt 12 is turbulent and the Navier-Stokes equation are connected and solved to solve the convection. A fifth step for further determining the internal temperature distribution of the melt 12 in consideration, a sixth step for determining the solid-liquid interface shape of the single crystal 14 and the melt 12 in accordance with an isotherm including the triple point S of the single crystal, Steps 4 to 6 are repeated until the triple point S reaches the melting point of the single crystal 14, the temperature distribution in the puller 11 is calculated, the mesh coordinates and the temperature of the single crystal 14 are obtained, and these data are respectively calculated. The seventh step input to the computer and the pulling length and height of the single crystal 14 are changed stepwise to calculate the temperature distribution in the pulling machine 11 repeatedly from the second step to the seventh step. Obtain the coordinates and temperature of the mesh and input these data to the computer, respectively, and the time from the start of pulling the single crystal 14 from the melt 12 to the completion of cooling of the single crystal 14 at predetermined intervals. Separation The ninth step of obtaining 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 obtained in the eighth step at each divided time interval. And vacancy and interstitial atoms after a 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. A tenth step for determining the concentration distribution of the first, an eleventh step for determining the formation start temperature of the void based on the concentration distribution of the vacancies, and a twelfth step for determining the formation start temperature of the high-temperature oxygen precipitate based on the concentration distribution of the vacancies. A step, a thirteenth step of obtaining a void concentration when the temperature at the lattice point of each mesh in the single crystal 14 is gradually decreased to a void formation start temperature, and a step of each mesh in the single crystal 14 A fourteenth step of determining the radius of the void when the temperature at the lattice point is lower than the void formation start temperature; and the temperature at the lattice point of each mesh in the single crystal 14 is gradually decreased to form a high temperature oxygen precipitate formation start temperature. A fifteenth step for determining the concentration of the high-temperature oxygen precipitates when the temperature of A sixteenth step for determining the radius of the high-temperature oxygen precipitate when it is lower than the formation start temperature of the elementary precipitate; a seventeenth step for repeating the ninth to sixteenth steps until the cooling of the single crystal 14 is completed; A first isoconcentration line HC1 that partitions a P band in which an object exists and a defect-free region in which neither a void nor a high-temperature oxygen precipitate exists and that the defect-free region is located on the single crystal bottom side from the P band_XAnd a first distribution line BC1 that partitions a B band in which high-concentration oxygen precipitates are present and a defect-free region, and that the defect-free region is positioned on the single crystal top side from the B band._XThe 19th step to obtain the first isoconcentration line HC1_XMaximum inflection point and first distribution line BC1_XThe difference ΔZ from the minimum value of the inflection point₁20th step for obtaining the parameter of the single crystal manufacturing condition P₂To P_NTo the first isoconcentration line HC1 by repeating the second step to the twentieth step._XMaximum inflection point and first distribution line BC1_XThe difference ΔZ from the minimum value of the inflection point₁After obtaining the difference ΔZ₁This is a simulation method for maximizing a defect-free region of a single crystal using a computer including a twenty-first step for obtaining a single crystal manufacturing condition that maximizes the single crystal.
[0010]
In the simulation method for maximizing the defect-free region of the single crystal according to claim 1, the parameter P₁In addition to obtaining the temperature distribution in the single crystal 14 grown from the melt 12 in consideration of the convection of the melt 12 under the single crystal production conditions, the temperature distribution in the single crystal 14 during the cooling process is also obtained. 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 concentration distribution and size distribution of voids in the single crystal 14 are predicted. In addition, the concentration distribution and size distribution of the high-temperature oxygen precipitates in the single crystal 14 are also predicted. Then, the first isoconcentration line HC1 in the single crystal 14_XAnd the first distribution line BC1_XIs obtained by calculation, and then the first isoconcentration line HC1_XMaximum inflection point and first distribution line BC1_XThe difference ΔZ from the minimum value of the inflection point₁Is calculated.
Next, the parameters of the above single crystal manufacturing conditions are set to P₂P in order from_NUntil the first isoconcentration line HC1_XMaximum inflection point and first distribution line BC1_XThe difference ΔZ from the minimum value of the inflection point₁Is calculated. Further, the first isoconcentration line HC1_XMaximum inflection point and first distribution line BC1_XThe difference ΔZ from the minimum value of the inflection point₁The condition for producing a single crystal with the largest value is obtained. Thereby, the single crystal manufacturing conditions under which the defect-free region is expanded most in the pulling direction and the radial direction of the single crystal 14 can be accurately obtained by calculation.
[0011]
In the present specification, the “high temperature oxygen precipitate” refers to a material in which the inner wall oxide film precipitated inside the void defect is treated like a sphere. The reason why the oxygen precipitate is limited to a high temperature is to distinguish it from a precipitate nucleus such as BMD (Bulk Micro Defect) that precipitates at a temperature lower than the void defect. In this specification, the “defect-free region” includes a defect-free region including an intrinsic point defect in which vacancies are dominant and a defect-free region including an intrinsic defect in which interstitial atoms are dominant.
[0012]
As shown in FIG. 6 and FIGS. 8 to 12, the invention according to claim 2 is a first step to a seventeenth step according to claim 1, a P band, a void and a high temperature oxygen in which high temperature oxygen precipitates are present. A second isoconcentration line HC2 that divides the defect-free region where none of the precipitates exist and divides the defect-free region so as to be positioned on the single crystal top side from the P band._XAnd a second distribution line BC2 that partitions the B band in which high-concentration oxygen precipitates are present from the defect-free region and that the defect-free region is positioned on the single crystal bottom side from the B band._XThe 19th step for obtaining the second isoconcentration line HC2_XMinimum inflection point and second distribution line BC2_XThe difference ΔZ from the maximum inflection point₂20th step for obtaining the parameter of the single crystal manufacturing condition P₂To P_NTo the second isoconcentration line HC2 by repeating the second step to the twentieth step._XMinimum inflection point and second distribution line BC2_XThe difference ΔZ from the maximum inflection point₂After obtaining the difference ΔZ₂This is a simulation method for maximizing a defect-free region of a single crystal using a computer including a twenty-first step for obtaining a single crystal manufacturing condition that maximizes the single crystal.
[0013]
In the simulation method for maximizing the defect-free region of the single crystal according to claim 2, the parameter P₁In the same manner as in the first aspect, the void concentration distribution and size distribution in the single crystal 14 are predicted, and the high-temperature oxygen precipitate concentration distribution and size distribution in the single crystal 14 are also predicted. Then, the second isoconcentration line HC2 in the single crystal 14_XAnd second distribution line BC2_XIs obtained by calculation, and then the second isoconcentration line HC2_XMaximum inflection point and second distribution line BC2_XThe difference ΔZ from the minimum value of the inflection point₂Is calculated.
Next, the parameters of the above single crystal manufacturing conditions are set to P₂P in order from_NUntil the second isoconcentration line HC2_XMaximum inflection point and second distribution line BC2_XThe difference ΔZ from the minimum value of the inflection point₂Is calculated. Further, the second isoconcentration line HC2_XMaximum inflection point and second distribution line BC2_XThe difference ΔZ from the minimum value of the inflection point₂The condition for producing a single crystal with the largest value is obtained. Thereby, the single crystal manufacturing conditions under which the defect-free region is expanded most in the pulling direction and the radial direction of the single crystal 14 can be accurately obtained by calculation.
[0014]
DETAILED DESCRIPTION OF THE INVENTION
Next, embodiments of the present invention will be described with reference to the drawings.
As shown in FIG. 6, a 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.
[0015]
A first simulation method for maximizing the defect-free region of the silicon single crystal 14 in the silicon single crystal puller 11 configured as described above will be described with reference to FIGS.
(1) Concentration distribution and size distribution of voids in the silicon single crystal 14
First, as a first step, single crystal manufacturing conditions for pulling the silicon single crystal 14 by the silicon single crystal puller 11 are arbitrarily defined. This single crystal manufacturing condition is a parameter group P in which a variable fed back to a hot zone of the pulling machine 11 described later is changed at regular intervals when the silicon single crystal 11 is pulled by the pulling machine 11.₁, P₂, ..., P_NIt is. The single crystal manufacturing conditions include the pulling speed of the silicon single crystal 14, the rotating speed of the silicon single crystal 14, the rotating speed of the quartz crucible 13, the flow rate of argon gas, the shape and material of the members constituting the heat cap, Examples include a gap between the lower end and the surface of the silicon melt 12, a heater output, and the like.
[0016]
As a second step, the silicon single crystal 14 has a predetermined length L₁Each member of the hot zone of the 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, and the like are modeled by mesh division. . 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.
[0017]
The reason why the radial mesh is limited to the range of 0.01 to 5.00 mm is that the calculation time becomes extremely long if it is less than 0.01 mm, and the calculation becomes unstable if it exceeds 5.00 mm. This is because the solid-liquid interface shape cannot be fixed. The longitudinal mesh is limited to the range of 0.01 to 5.00 mm because the calculation time is extremely long if it is less than 0.01 mm, and the calculated value of the solid-liquid interface shape matches the actual measurement value if it exceeds 5.00 mm. Because it will not do. When a part of the radial mesh is limited to the range of 0.01 to 5.00, the 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.
[0018]
As a third step, the meshes are grouped for each member in the hot zone, and the physical property values of the members are input to the computer for the grouped meshes. For example, if the chamber is made of stainless steel, the thermal conductivity, emissivity, viscosity, volume expansion coefficient, density and specific heat of the stainless steel are input to the computer. Further, the pulling length of the silicon 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.
[0019]
As a fourth step, the surface temperature distribution of each member in the hot zone is obtained by calculation using a computer based on the amount of heat generated by the heater and the radiation rate of each member. That is, the heating value of the heater is arbitrarily set and input to the computer, and the surface temperature distribution of each member is obtained by calculation using the computer from the emissivity of each member. Next, as a fifth step, the internal temperature distribution of each member is calculated by solving the heat conduction equation (1) using a computer based on the surface temperature distribution and the thermal conductivity of each member in the hot zone. Here, in order to simplify the description, an xyz orthogonal coordinate system is used, but a cylindrical coordinate system is used in actual calculation.
[0020]
[Expression 1]

Where ρ 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, 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.
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.
[0021]
[Expression 2]

Where κ_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.
[0022]
[Equation 3]

Here, 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.
[0023]
The turbulent model equation (2) is referred to as a kl-model equation, and an arbitrary value within the range of 0.4 to 0.6 is preferably used as the turbulent parameter C of the model equation. The reason why the turbulent flow parameter C is limited to the range of 0.4 to 0.6 is that when the value is less than 0.4 or exceeds 0.6, there is a problem that the interface shape obtained by calculation does not match the actual measurement value. . The Navier-Stokes equations (3) to (5) are equations of motion when the silicon melt 12 is a fluid that is incompressible and has a constant viscosity.
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. .
[0024]
[Expression 4]

Here, u, v, and w are the x, y, and z direction components of the flow velocity at each mesh point of the silicon melt 12, T is the absolute temperature at each mesh point of the silicon melt 12, and ρ is The density of the silicon melt 12, c is the specific heat of the silicon melt 12, and κ_lIs the molecular thermal conductivity (physical property value) and κ_tIs the turbulent thermal conductivity calculated using equation (1).
[0025]
Next, as a sixth step, the silicon-crystalline interface 14 and the silicon melt 12 include a silicon triple point S (solid-liquid-gas triple point) indicated by a point S in FIG. It is obtained by calculation using a computer according to the isotherm. As a seventh step, the amount of heat generated by the heater input to the computer is changed (increase gradually), and the fourth to sixth steps are repeated until the triple point S reaches the melting point of the silicon single crystal 14, and then the puller The temperature distribution in 11 is calculated to obtain the coordinates and temperature of the mesh of the silicon single crystal, and these data are stored in the computer.
[0026]
Next, as an eighth step, the pulling length L of the silicon single crystal 14₁After adding δ (for example, 50 mm) to the above and repeating the second to seventh steps, 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 the eighth 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. 6). ) 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.₁And δ (for example, 50 mm) is added, and the second to seventh steps are repeated as described above.
[0027]
Lifting height H of silicon single crystal 14₁Is H₂When the value is reached, the process proceeds to the ninth step. In the ninth 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 constants used in equations for determining the concentration distribution and size distribution of voids and high-temperature oxygen precipitates described later are respectively calculated by computers. To enter. 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 eighth step at every divided time interval Δt seconds.₁And pulling height H₁And the temperature distribution in the silicon single crystal 14 are calculated.
[0028]
That is, in the second to eighth steps, the coordinates and temperature of the silicon single crystal mesh are obtained for each pulling length δ, and it takes several tens of minutes to pull the silicon single crystal, for example, 50 mm. By differentiating the temperature change of the single crystal mesh as a function of time, the time 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, 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, the concentration distribution of the vacancies and interstitial silicon after lapse of Δt seconds is obtained by calculation (No. 1). 10 steps).
[0029]
Specifically, the vacancy concentration C_vIs calculated by the following equation (7), and the interstitial silicon concentration C_iIs calculated by the following equation (8). In formula (7) and formula (8), concentration C_vAnd concentration 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.
[0030]
[Equation 5]

Where K₁And K₂Is a constant, E_iAnd E_vAre the energy of formation of interstitial silicon and vacancies, respectively, and C_v ^eAnd C_i ^eThe letter “e” of “” means the equilibrium amount, “k” means Boltzmann constant, and “T” means absolute temperature.
The above equilibrium equation is differentiated with respect to time and becomes the following equations (9) and (10) for the vacancies and interstitial silicon, respectively.
[0031]
[Formula 6]

Here, Θ (x) is a heaviside function. That is, Θ (x) = 0 when x <0, and Θ (x) = 1 when x> 0. T_nIs the formation temperature T of the high temperature oxygen precipitate_pAnd void formation start temperature T_v0Is the higher temperature when compared with. 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 represented by the following equations (11) and (12).
[0032]
[Expression 7]

Where △ E_vAnd △ E_iAre the activation energies of the vacancies and interstitial silicon, respectively, d_vAnd d_iAre constants. Further, in the second term on the right side of each of the equations (9) and (10),
[Equation 8]

Is 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)_v0Is the void concentration, r_v0Is the radius of the void, and N in the fifth term on the right side of equation (9)_pIs the concentration of high temperature oxygen precipitates and R_pIs the radius of the high temperature oxygen precipitate and γ is SiO₂This is the vacancy consumption amount per oxygen atom required for the precipitate to be deposited without distortion.
[0033]
On the other hand, the above formula (9) is satisfied because the flow velocity for the precipitation of the vacancies is sufficiently large, and the Si matrix and SiO₂And the volume difference per unit mass can be filled, that is, D_v(C_v-C_v ^e) ≧ γD₀C₀This is the case. In cases other than the above, the following equation (13) is established.
[0034]
[Equation 9]

Next, as an eleventh step, the vacancy concentration C determined by solving the above diffusion equation_vBased on the distribution, void formation start temperature T_v0Is obtained from the following equation (14).
[0035]
[Expression 10]

Where C_vmIs the melting point T of the silicon melt_mIs the vacancy equilibrium concentration at E_vIs the vacancy formation energy and T_mIs the melting point temperature of the silicon single crystal 14. 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.
[0036]
The twelfth step will be described later, and the thirteenth step will be described. In the thirteenth step, the temperature at the lattice point of each mesh in the silicon single crystal 14 gradually decreases, and the void formation start temperature T_v0The void concentration N using the following approximate expression (15):_v0Ask for.
[0037]
[Expression 11]

Further, as a fourteenth step, the temperature at the lattice point of each mesh in the silicon single crystal 14 is the void formation start temperature T._v0Lower void radius r_v0Is obtained from the following equation (16).
[0038]
[Expression 12]

Where t₁The temperature at the lattice point of the mesh of the silicon single crystal 14 is the void formation start temperature T_v0Is the time when it drops to r_crIs the critical diameter of the void. The ninth to eleventh steps, the thirteenth and the fourteenth steps are repeated until the silicon single crystal 14 is cooled to a specific value between 800 and 1000 ° C., for example, 900 ° C. or less (17th step). Equations (9) to (16) are coupled and solved by a computer. The fifteenth and sixteenth steps will be described later.
[0039]
(2) Concentration distribution and size distribution of high-temperature oxygen precipitates in the silicon single crystal 14
Returning to the twelfth step, the vacancy concentration C determined by solving the diffusion equations (9) and (10)._vBased on the distribution, the formation start temperature T of the high-temperature oxygen precipitate_pIs calculated from the following equation (17).
[0040]
[Formula 13]

Where C₀Is the oxygen concentration and C_0mIs the melting point T of the silicon melt 12_mIs the oxygen equilibrium concentration at E₀Is oxygen dissolution energy. E_vIs the vacancy formation energy and σ_pIs Si and SiO in the silicon single crystal 14₂And the interfacial energy. Further, γ is SiO₂This is the vacancy consumption amount per oxygen atom necessary for the precipitate to be deposited without distortion, and its value is 0.68.
[0041]
Next, as a fifteenth 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 high-temperature oxygen precipitates_pThe concentration N of the high-temperature oxygen precipitate is calculated using the following approximate equation (18)._pAsk for. In formula (18), a₂Is a constant.
[0042]
[Expression 14]

Furthermore, as a sixteenth step, the temperature at the lattice point of each mesh in the silicon single crystal 14 is the temperature at which the high temperature oxygen precipitate formation start T_pRadius R of hot oxygen precipitate at lower_pIs obtained from the following equation (19).
[0043]
[Expression 15]

Where t₂The temperature at the lattice point of the mesh of the silicon single crystal 14 is the formation start temperature T of the high temperature oxygen precipitate._pIs the time when it drops to R,_crIs the critical diameter of the high temperature oxygen precipitate.
[0044]
On the other hand, the above equation (19) is satisfied because the flow rate for the precipitation of vacancies is sufficiently large.₂And the volume difference per unit mass can be filled, that is, D_v(C_v-C_v ^e) ≧ γD₀C₀This is the case. In cases other than the above, the following equation (20) is established.
[0045]
[Expression 16]

The ninth, tenth, twelfth, fifteenth and sixteenth steps are repeated until the silicon single crystal 14 cools to a specific value between 800 and 1000 ° C., for example, 900 ° C. or less (17th step). The above equations (9) to (13) and equations (17) to (20) are coupled and solved by a computer.
[0046]
As an eighteenth step, the P band where the high-temperature oxygen precipitates are present and a defect-free region where neither voids nor high-temperature oxygen precipitates exist are defined, and the defect-free region is located on the single crystal bottom side from the P band. First concentration line HC1_XIs obtained by calculation (FIG. 7). Next, as a nineteenth step, a first distribution line BC1 that partitions a B band in which high-concentration oxygen precipitates are present from a defect-free region and that the defect-free region is positioned on the single crystal top side from the B band._XIs calculated. The first distribution line is an isoconcentration line in which the interstitial silicon concentration is a specific value in the range of 0.12 to 0.13% relative to the equilibrium concentration of interstitial silicon at the silicon melting point, for example, 0.126%. Say. Here, the specific value in the range of the interstitial silicon concentration of 0.12 to 0.13%, for example, 0.126%, is observed from the lifetime map or the like among the interstitial silicon concentrations obtained by performing the simulation. The density corresponds to the B band. The B band refers to a region in which an aggregate of interstitial silicon serves as a nucleus and oxygen precipitation occurs at a high concentration by heat treatment.
[0047]
Thereafter, as the 20th step, the first isoconcentration line HC1_XMaximum inflection point and first distribution line BC1_XThe difference ΔZ from the minimum value of the inflection point₁Find (X). First isoconcentration line HC1_XIf there are, for example, three inflection points, the coordinates of each inflection point Q1 (1, X), Q1 (2, X) and Q1 (3, X) are expressed as (r_Q1(1, X), z_Q1(1, X)), (r_Q1(2, X), z_Q1(2, X)) and (r_Q1(3, X), z_Q1(3, X)) of these coordinates, z coordinate z_Q1(1, X), z_Q1(2, X) and z_Q1The maximum value z of (3, X)_Q1Find (1, X). The first distribution line BC1_XFor example, when there are five inflection points, each inflection point S1 (1, X), S1 (2, X), S1 (3, X), S1 (4, X) and S1 (5, X) ) Coordinate (r_S1(1, X), z_S1(1, X)), (r_S1(2, X), z_S1(2, X)), (r_S1(3, X), z_S1(3, X)), (r_S1(4, X), z_S1(4, X)) and (r_S1(5, X), z_S1(5, X)) of these coordinates, z coordinate z_S1(1, X), z_S1(2, X), z_S1(3, X), z_S1(4, X) and z_S1Minimum value z of (5, X)_S1Find (3, X). And the first isoconcentration line HC1_XMaximum value of inflection point z_Q1(1, X) and first distribution line BC1_XMinimum value z of inflection point_S1ΔZ difference from (3, X)₁Find (X).
[0048]
As the 21st step, the parameters of the single crystal manufacturing conditions are set to P₂Execute steps 2 through 20 to change the parameters of the single crystal manufacturing conditions to P_ThreeThe parameters of the single crystal manufacturing conditions are P, such as executing the second step to the twentieth step instead of_NAfter repeating the second step to the twentieth step until the first concentration line HC1 is reached_XMaximum inflection point and first distribution line BC1_XThe difference ΔZ from the minimum value of the inflection point₁(X) is calculated, and this difference ΔZ₁A single crystal production condition in which (X) is the largest is obtained by calculation.
[0049]
The single crystal manufacturing conditions in the silicon single crystal pulling machine 11 calculated in this way are the conditions that the defect-free region expands most in the pulling direction and the radial direction of the silicon single crystal 14. It can be accurately obtained by the above calculation. As a result, it is possible to set a single crystal manufacturing condition in which the defect-free region of the silicon single crystal 14 is expanded most by calculation using a computer as described above, so the lifetime is measured by actually pulling up the silicon single crystal. The experiment to be performed is only necessary to confirm whether or not the defect-free region is actually expanded under the above-mentioned single crystal manufacturing conditions. Therefore, since the number of experiments can be reduced, optimal single crystal production conditions can be obtained with relatively little time and effort. As is clear from FIG. 7, the difference ΔZ₁As (X) increases, the width of the pulling speed of the silicon single crystal 14 for forming a defect-free region in the pulling direction and the radial direction of the silicon single crystal 14, that is, the pure margin increases, so that there is no defect over the entire length. The silicon single crystal 14 can be pulled up relatively easily.
[0050]
Next, a second simulation method for maximizing the defect-free region of the silicon single crystal 14 in the silicon single crystal puller 11 will be described with reference to FIGS.
The first to seventeenth steps are the same as the first to seventeenth steps of the first simulation method.
[0051]
As an eighteenth step, the P band where the high-temperature oxygen precipitates are present and the defect-free region where neither voids nor high-temperature oxygen precipitates exist are defined, and the defect-free region is located on the single crystal top side from the P band. Second isoconcentration line HC2_XIs obtained by calculation (FIG. 7). Next, as a nineteenth step, a second distribution line BC2 that partitions the B band in which high-concentration oxygen precipitates are present and the defect-free region, and partitions the defect-free region so that it is located on the single crystal bottom side from the B band._XIs calculated. The second distribution line is an isoconcentration line in which the interstitial silicon concentration is a specific value in the range of 0.12 to 0.13% relative to the equilibrium concentration of interstitial silicon at the silicon melting point, for example, 0.126%. Say. Here, the specific value in the range of the interstitial silicon concentration of 0.12 to 0.13%, for example, 0.126%, is observed from the lifetime map or the like among the interstitial silicon concentrations obtained by performing the simulation. The density corresponds to the B band.
[0052]
Thereafter, as the 20th step, the second isoconcentration line HC2_XMinimum inflection point and second distribution line BC2_XThe difference ΔZ from the maximum inflection point₂Find (X). Second isoconcentration line HC2_XIf there are, for example, three inflection points, the coordinates of each inflection point Q2 (1, X), Q2 (2, X) and Q2 (3, X) are expressed as (r_Q2(1, X), z_Q2(1, X)), (r_Q2(2, X), z_Q2(2, X)) and (r_Q2(3, X), z_Q2(3, X)) of these coordinates, z coordinate z_Q2(1, X), z_Q2(2, X) and z_Q2Minimum value z of (3, X)_Q2Find (1, X). The second distribution line BC2_XFor example, when there are five inflection points, each inflection point S2 (1, X), S2 (2, X), S2 (3, X), S2 (4, X) and S2 (5, X ) Coordinate (r_S2(1, X), z_S2(1, X)), (r_S2(2, X), z_S2(2, X)), (r_S2(3, X), z_S2(3, X)), (r_S2(4, X), z_S2(4, X)) and (r_S2(5, X), z_S2(5, X)) of these coordinates, z coordinate z_S2(1, X), z_S2(2, X), z_S2(3, X), z_S2(4, X) and z_S2The maximum value z of (5, X)_S2Find (3, X). And the second isoconcentration line HC2_XMinimum value z of inflection point_Q2(1, X) and second distribution line BC2_XMaximum value of inflection point z_S2ΔZ difference from (3, X)₂Find (X).
[0053]
As the 21st step, the parameters of the single crystal manufacturing conditions are set to P₂Execute steps 2 through 20 to change the parameters of the single crystal manufacturing conditions to P_ThreeThe parameters of the single crystal manufacturing conditions are P, such as executing the second step to the twentieth step instead of_NAfter repeating the second step to the twentieth step until the second concentration line HC2_XMinimum inflection point and second distribution line BC2_XThe difference ΔZ from the maximum inflection point₂(X) is calculated, and this difference ΔZ₂A single crystal production condition in which (X) is the largest is obtained by calculation.
[0054]
The single crystal manufacturing conditions in the silicon single crystal pulling machine 11 calculated in this way are the conditions that the defect-free region expands most in the pulling direction and the radial direction of the silicon single crystal 14. It can be accurately obtained by the above calculation. As a result, it is possible to set a single crystal manufacturing condition in which the defect-free region of the silicon single crystal 14 is expanded most by calculation using a computer as described above, so the lifetime is measured by actually pulling up the silicon single crystal. The experiment to be performed is only necessary to confirm whether or not the defect-free region is actually expanded under the above-mentioned single crystal manufacturing conditions. Therefore, since the number of experiments can be reduced, optimal single crystal production conditions can be obtained with relatively little time and effort. As is clear from FIG. 7, the difference ΔZ₂As (X) increases, the width of the pulling speed of the silicon single crystal 14 for forming a defect-free region in the pulling direction and the radial direction of the silicon single crystal 14, that is, the pure margin increases, so that there is no defect over the entire length. The silicon single crystal 14 can be pulled up relatively easily.
In the above embodiment, a silicon single crystal is used as the single crystal. However, a GaAs single crystal, an InP single crystal, a ZnS single crystal, or a ZnSe single crystal may be used.
[0055]
【The invention's effect】
As described above, according to the present invention, by first considering the convection of the melt under predetermined single crystal manufacturing conditions, the calculated crystal and the solid-liquid interface shape of the melt can be obtained when the actual single crystal is pulled up. In addition to determining the internal temperature of the single crystal so that it substantially matches the shape, it is also possible to determine the temperature distribution within the single crystal during the cooling process, that is, the single crystal in the cooling process of the single crystal separated from the melt. By considering the effect of slow cooling and rapid cooling of the crystal, the concentration distribution and size distribution of defects in the single crystal are obtained using a computer. Next, the first isoconcentration line and the first distribution line in the single crystal are obtained by calculation, and the difference between the maximum value of the inflection point of the first isoconcentration line and the minimum value of the inflection point of the first distribution line is calculated. After obtaining by calculation, the parameters of the single crystal production conditions are changed, and the difference between the maximum value of the inflection point of the first isoconcentration line and the minimum value of the inflection point of the first distribution line is calculated in the same manner as described above. Ask for. Further, a single crystal production condition is obtained in which the difference between the maximum value of the inflection point of the first isoconcentration line and the minimum value of the inflection point of the first distribution line is the largest. As a result, the single crystal manufacturing conditions under which the defect-free region expands most in the pulling direction and the radial direction of the single crystal can be accurately obtained by calculation. For this reason, compared with a conventional defect-free single crystal ingot that required many experiments to set a single crystal production condition in which the defect-free region is most expanded, in the present invention, the defect-free region is The number of experiments for setting the most expanded single crystal production conditions can be reduced. As a result, in the present invention, it is possible to obtain a single crystal manufacturing condition in which the defect-free region is most expanded with relatively little time and effort.
[0056]
Also, the second isoconcentration line and the second distribution line in the single crystal are obtained by calculation, and the difference between the maximum value of the inflection point of the second isoconcentration line and the minimum value of the inflection point of the second distribution line is calculated. Then, the parameter of the single crystal production condition is changed, and the difference between the maximum value of the inflection point of the second isoconcentration line and the minimum value of the inflection point of the second distribution line is calculated in the same manner as above. Further, the same effect as described above can be obtained by obtaining a single crystal production condition in which the difference between the maximum value of the inflection point of the first isoconcentration line and the minimum value of the inflection point of the first distribution line is maximized. can get.
[Brief description of the drawings]
FIG. 1 is a flowchart showing a first stage of a first simulation method according to an embodiment of the present invention.
FIG. 2 is a flowchart showing a second stage of the simulation method.
FIG. 3 is a flowchart showing a third stage of the simulation method.
FIG. 4 is a flowchart showing a fourth stage of the simulation method.
FIG. 5 is a flowchart showing a fifth stage of the simulation method.
FIG. 6 is a cross-sectional view of an essential part of a silicon single crystal pulling machine having a silicon melt mesh structure.
FIG. 7 is an explanatory diagram showing the distribution of interstitial silicon and vacancies in the silicon single crystal when the pulling rate of the silicon single crystal is changed.
FIG. 8 is a flowchart showing a first stage of a second simulation method according to the embodiment of the present invention.
FIG. 9 is a flowchart showing a second stage of the simulation method.
FIG. 10 is a flowchart showing a third stage of the simulation method.
FIG. 11 is a flowchart showing a fourth stage of the simulation method.
FIG. 12 is a flowchart showing the fifth stage of the simulation method.
[Explanation of symbols]
11 Silicon single crystal pulling machine
12 Silicon melt
14 Silicon single crystal
S Triple point of silicon

Claims (2)

Parameter groups P ₁ , P ₂ ,..., P _N are arbitrarily defined by varying variables fed back as a single crystal production condition when pulling up the single crystal 14 by the puller 11. The first step;
From the parameter group P ₁ , P ₂ ,..., P _N , the parameter P ₁ is used as the single crystal production condition, and the cooling of the single crystal (14) is completed from the start of the pulling of the single crystal (14) from the melt. A second step of modeling the hot zone of the puller (11) up to a mesh structure;
The summarized mesh for each member of the hot zone and the with physical properties of each member for the gathered mesh parameter groups P _1, P _2, ..., first enter the parameters P ₁ of the P _N to the computer 3 steps,
A fourth step of obtaining the surface temperature distribution of each member based on the heat value of the heater and the radiation rate of each member;
The turbulence obtained by assuming that the melt (12) is turbulent after obtaining the internal temperature distribution of each member by solving the heat conduction equation based on the surface temperature distribution and thermal conductivity of each member. A fifth step of further determining the internal temperature distribution of the melt (12) considering convection by connecting and solving the flow model equation and the Navier-Stokes equations;
A sixth step of obtaining a solid-liquid interface shape of the single crystal (14) and the melt (12) according to an isotherm including a triple point (S) of the single crystal;
The fourth step to the sixth step are repeated until the triple point (S) reaches the melting point of the single crystal (14), and the temperature distribution in the puller (11) is calculated to calculate the single crystal (14). A seventh step of obtaining the coordinates and temperature of the mesh and inputting these data to the computer, respectively;
The pulling length and the pulling height of the single crystal (14) are changed stepwise, and the temperature distribution in the pulling machine (11) is calculated repeatedly from the second step to the seventh step to calculate the single crystal (14 8) to obtain the coordinates and temperature of the mesh and input these data to the computer,
Eighth step of dividing the time from the start of pulling of the single crystal (14) from the melt (12) to the completion of cooling of the single crystal (14) at predetermined intervals. A ninth step of determining the pulling length and height of the single crystal (14) and the temperature distribution in the single crystal (14) from the mesh coordinates and temperature data of the single crystal (14) determined in
The concentration distribution of vacancies and interstitial atoms after the predetermined time interval has been obtained by solving the diffusion equation based on the diffusion coefficients and boundary conditions of vacancies and interstitial atoms in the single crystal (14) 10th step;
An eleventh step of determining a void formation start temperature based on the concentration distribution of the pores;
A twelfth step of determining a formation start temperature of the high-temperature oxygen precipitate based on the concentration distribution of the pores;
A thirteenth step of determining the concentration of the void when the temperature at the lattice point of each mesh in the single crystal (14) gradually decreases to the formation start temperature of the void;
A fourteenth step of determining a radius of the void when a temperature at a lattice point of each mesh in the single crystal (14) is lower than a formation start temperature of the void;
A fifteenth step of determining the concentration of the high-temperature oxygen precipitates when the temperature at the lattice points of each mesh in the single crystal (14) gradually decreases to reach the formation start temperature of the high-temperature oxygen precipitates;
A sixteenth step of obtaining a radius of the high-temperature oxygen precipitate when the temperature at the lattice point of each mesh in the single crystal (14) is lower than the formation start temperature of the high-temperature oxygen precipitate;
A seventeenth step of repeating the ninth to sixteenth steps until the cooling of the single crystal (14) is completed;
The P band where the high-temperature oxygen precipitates are present and the defect-free region where neither the void nor the high-temperature oxygen precipitates exist are partitioned, and the defect-free region is located on the single crystal bottom side from the P band. An eighteenth step for determining a first isoconcentration line (HC1 _X ) to be partitioned;
A first distribution line (BC1 _X ) is obtained that partitions the B band in which high-concentration oxygen precipitates are present and the defect-free region, and partitions the defect-free region so that it is located on the single crystal top side from the B band. 19th step;
A twentieth step of obtaining a difference (ΔZ ₁ ) between the maximum value of the inflection point of the first isoconcentration line (HC1 _X ) and the minimum value of the inflection point of the first distribution line (BC1 _X );
The parameter of the single crystal production condition is sequentially changed from P ₂ to _PN, and the second step to the twentieth step are repeated, so that the maximum value of the inflection point of the first isoconcentration line (HC1 _X ) and the first distribution are obtained. A computer using a computer including a twenty-first step of obtaining a single crystal production condition in which the difference (ΔZ ₁ ) is largest after obtaining a difference (ΔZ ₁ ) from the minimum value of the inflection point of the line (BC1 _X ). A simulation method that maximizes the defect-free region of a crystal. Parameter groups P ₁ , P ₂ ,..., P _N are arbitrarily defined by varying variables fed back as a single crystal production condition when pulling up the single crystal 14 by the puller 11. The first step;
From the parameter group P ₁ , P ₂ ,..., P _N , the parameter P ₁ is used as the single crystal production condition, and the cooling of the single crystal (14) is completed from the start of the pulling of the single crystal (14) from the melt. A second step of modeling the hot zone of the puller (11) up to a mesh structure;
The summarized mesh for each member of the hot zone and the with physical properties of each member for the gathered mesh parameter groups P _1, P _2, ..., first enter the parameters P ₁ of the P _N to the computer 3 steps,
A fourth step of obtaining the surface temperature distribution of each member based on the heat value of the heater and the radiation rate of each member;
The turbulence obtained by assuming that the melt (12) is turbulent after obtaining the internal temperature distribution of each member by solving the heat conduction equation based on the surface temperature distribution and thermal conductivity of each member. A fifth step of further determining the internal temperature distribution of the melt (12) considering convection by connecting and solving the flow model equation and the Navier-Stokes equations;
A sixth step of obtaining a solid-liquid interface shape of the single crystal (14) and the melt (12) according to an isotherm including a triple point (S) of the single crystal;
The fourth step to the sixth step are repeated until the triple point (S) reaches the melting point of the single crystal (14), and the temperature distribution in the puller (11) is calculated to calculate the single crystal (14). A seventh step of obtaining the coordinates and temperature of the mesh and inputting these data to the computer, respectively;
The pulling length and the pulling height of the single crystal (14) are changed stepwise, and the temperature distribution in the pulling machine (11) is calculated repeatedly from the second step to the seventh step to calculate the single crystal (14 8) to obtain the coordinates and temperature of the mesh and input these data to the computer,
Eighth step of dividing the time from the start of pulling of the single crystal (14) from the melt (12) to the completion of cooling of the single crystal (14) at predetermined intervals. A ninth step of determining the pulling length and height of the single crystal (14) and the temperature distribution in the single crystal (14) from the mesh coordinates and temperature data of the single crystal (14) determined in
The concentration distribution of vacancies and interstitial atoms after the predetermined time interval has been obtained by solving the diffusion equation based on the diffusion coefficients and boundary conditions of vacancies and interstitial atoms in the single crystal (14) 10th step;
An eleventh step of determining a void formation start temperature based on the concentration distribution of the pores;
A twelfth step of determining a formation start temperature of the high-temperature oxygen precipitate based on the concentration distribution of the pores;
A thirteenth step of determining the concentration of the void when the temperature at the lattice point of each mesh in the single crystal (14) gradually decreases to the formation start temperature of the void;
A fourteenth step of determining a radius of the void when a temperature at a lattice point of each mesh in the single crystal (14) is lower than a formation start temperature of the void;
A fifteenth step of determining the concentration of the high-temperature oxygen precipitates when the temperature at the lattice points of each mesh in the single crystal (14) gradually decreases to reach the formation start temperature of the high-temperature oxygen precipitates;
A sixteenth step of obtaining a radius of the high-temperature oxygen precipitate when the temperature at the lattice point of each mesh in the single crystal (14) is lower than the formation start temperature of the high-temperature oxygen precipitate;
A seventeenth step of repeating the ninth to sixteenth steps until the cooling of the single crystal (14) is completed;
The P band where the high-temperature oxygen precipitates are present and the defect-free region where neither the void nor the high-temperature oxygen precipitates exist are partitioned, and the defect-free region is located on the single crystal top side from the P band. An 18th step for obtaining a second isoconcentration line (HC2 _X ) to be partitioned;
Obtain a second distribution line (BC2 _X ) that partitions the B band in which high-concentration oxygen precipitates are present and the defect-free region and that the defect-free region is positioned on the single crystal bottom side from the B band. 19th step;
A twentieth step of obtaining a difference (ΔZ ₂ ) between the minimum value of the inflection point of the second isoconcentration line (HC2 _X ) and the maximum value of the inflection point of the second distribution line (BC2 _X );
The parameter of the single crystal production condition is sequentially changed from P ₂ to _PN, and the second step to the twentieth step are repeated, so that the minimum value of the inflection point of the second isoconcentration line (HC2 _X ) and the second distribution are obtained. A computer including a twenty-first step for obtaining a single crystal manufacturing condition in which the difference (ΔZ ₂ ) is the largest after obtaining the difference (ΔZ ₂ ) from the maximum value of the inflection point of the line (BC2 _X ). A simulation method that maximizes the defect-free region of a crystal.

Images