JP4154936B2

JP4154936B2 - Single crystal defect-free region simulation method

Info

Publication number: JP4154936B2
Application number: JP2002185046A
Authority: JP
Inventors: 浩之介北村; 純古川; 直樹小野
Original assignee: Sumco Corp
Current assignee: Sumco Corp
Priority date: 2002-06-25
Filing date: 2002-06-25
Publication date: 2008-09-24
Anticipated expiration: 2022-06-25
Also published as: JP2004026567A

Description

【０００１】
【発明の属する技術分野】
本発明は、チョクラルスキー（以下、ＣＺという。）法にて引上げられるシリコン等の単結晶の無欠陥領域をコンピュータシミュレーションする方法に関するものである。
【０００２】
【従来の技術】
近年の半導体集積回路の超微細化にともないデバイスの歩留まりを低下させる要因として、結晶に起因したパーティクル（Crystal Originated Particle、以下、ＣＯＰという。）や、酸化誘起積層欠陥（Oxidation Induced Stacking Fault、以下、ＯＩＳＦという。）の核となる酸素析出物の微小欠陥や、或いは侵入型転位（Interstitial-type Large Dislocation、以下、Ｌ／Ｄという。）の存在が挙げられている。
【０００３】
ＣＯＰは、鏡面研磨されたシリコンウェーハをアンモニアと過酸化水素の混合液で洗浄すると、ウェーハ表面に出現するピットである。このウェーハをパーティクルカウンタで測定すると、このピットがパーティクルとして検出される。このピットは結晶に起因する。ＣＯＰは電気的特性、例えば酸化膜の経時絶縁破壊特性（Time Dependent dielectric Breakdown、ＴＤＤＢ）、酸化膜耐圧特性（Time Zero Dielectric Breakdown、ＴＺＤＢ）等を劣化させる原因となる。またＣＯＰがウェーハ表面に存在するとデバイスの配線工程において段差を生じ、断線の原因となり得る。そして素子分離部分においてもリーク等の原因となり、製品の歩留りを低くする。
【０００４】
ＯＩＳＦは、結晶成長時に微小な酸素析出核がシリコン単結晶中に導入され、半導体デバイスを製造する際の熱酸化工程等で顕在化する欠陥である。このＯＩＳＦは、デバイスのリーク電流を増加させる等の不良原因になる。Ｌ／Ｄは、転位クラスタとも呼ばれたり、或いはこの欠陥を生じたシリコンウェーハをフッ酸を主成分とする選択エッチング液に浸漬するとピットを生じることから転位ピットとも呼ばれる。このＬ／Ｄも、電気的特性、例えばリーク特性、アイソレーション特性等を劣化させる原因となる。以上のことから、半導体集積回路を製造するために用いられるシリコンウェーハからＣＯＰ、ＯＩＳＦ及びＬ／Ｄを減少させることが必要となっている。
【０００５】
このＣＯＰ、ＯＩＳＦ及びＬ／Ｄを有しない無欠陥の単結晶インゴット及びこの単結晶インゴットからスライスされたシリコンウェーハが米国特許番号６，０４５，６１０号に対応する特開平１１−１３９３号公報に開示されている。この無欠陥の単結晶インゴットは、単結晶インゴット内での空孔型点欠陥の凝集体及び格子間シリコン型点欠陥の凝集体がそれぞれ存在しないパーフェクト領域を［Ｐ］とするとき、パーフェクト領域［Ｐ］からなる単結晶インゴットである。パーフェクト領域［Ｐ］は、単結晶インゴット内で空孔型点欠陥が優勢であって凝集した点欠陥を有する領域［Ｖ］と、格子間シリコン型点欠陥が優勢であって凝集した点欠陥を有する領域［Ｉ］との間に介在する。
【０００６】
【発明が解決しようとする課題】
しかし、上記従来の特開平１１−１３９３号公報に示された無欠陥の単結晶インゴットでは、単結晶インゴットの全長にわたって空孔型点欠陥の凝集体及び格子間シリコン型点欠陥の凝集体が存在しない引上げ条件を設定するために、引上げ機の種類や単結晶インゴットの太さ毎に実際に単結晶インゴットを引上げてそのライフタイムを測定するという実験を数多く繰返さなければならず、最適な引上げ条件を設定するのに極めて多くの時間と労力を要する不具合があった。
また、上記従来の無欠陥の単結晶インゴットでは、点欠陥濃度を求めるために用いているＶ／Ｇ理論、即ち単結晶の引上げ速度をＶ（ｍｍ／分）とし、単結晶及び融液の固液界面近傍の温度勾配をＧ（℃／ｍｍ）とするときに、Ｖ／Ｇ（ｍｍ²／分・℃）を制御するボロンコフの理論を簡略化した一次元モデル（単結晶の引上げ方向だけのモデル）であるため、単結晶の半径方向への点欠陥の拡散（外方拡散）が考慮されておらず、単結晶の半径方向への点欠陥の濃度勾配を求めることができない問題点があった。
更に、上記従来の無欠陥の単結晶インゴットでは、Ｖ／Ｇ理論にて点欠陥の濃度を求める場合に、上記Ｇを高精度で求めるために融液の対流をモデル化して計算する必要があるにも拘らず、このような考慮が全くなされていないため、単結晶内の点欠陥の濃度分布を精度良く求めることができない問題点もあった。
本発明の目的は、実際に単結晶を引上げるという実験の回数を少なくすることができ、これにより比較的少ない時間と労力で、単結晶の最適な引上げ条件を設定できる、単結晶の無欠陥領域シミュレーション方法を提供することにある。
【０００７】
【課題を解決するための手段】
請求項１に係る発明は、図１〜図５に示すように、引上げ機１１により単結晶１４を引上げるときの単結晶製造条件としてフィードバックされる変数を一定の間隔で変量させたパラメータ群Ｐ₁，Ｐ₂，…，Ｐ_Nを任意に定義する第１ステップと、引上げ機１１により単結晶１４を所定長さまで引上げた状態における単結晶１４の引上げ機１１のホットゾーンをメッシュ構造でモデル化する第２ステップと、ホットゾーンの各部材毎にメッシュをまとめかつこのまとめられたメッシュに対する各部材の物性値とともにパラメータ群Ｐ₁，Ｐ₂，…，Ｐ_NのうちのパラメータＰ₁をコンピュータに入力する第３ステップと、各部材の表面温度分布をヒータの発熱量及び各部材の輻射率に基づいて求める第４ステップと、各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式を解くことにより各部材の内部温度分布を求めた後に融液１２が乱流であると仮定して得られた乱流モデル式及びナビエ・ストークスの方程式を連結して解くことにより対流を考慮した融液１２の内部温度分布を更に求める第５ステップと、単結晶１４及び融液１２の固液界面形状を単結晶の三重点Ｓを含む等温線に合せて求める第６ステップと、第４ステップから第６ステップを三重点Ｓが単結晶１４の融点になるまで繰返し引上げ機１１内の温度分布を計算して単結晶１４のメッシュの座標及び温度を求めこれらのデータをそれぞれコンピュータに入力する第７ステップと、単結晶１４の引上げ長を段階的に変えて第２ステップから第７ステップまでを繰返し引上げ機１１内の温度分布を計算して単結晶１４のメッシュの座標及び温度を求めこれらのデータをそれぞれコンピュータに入力する第８ステップと、単結晶１４のメッシュの座標及び温度のデータと単結晶１４内の格子間原子及び空孔の拡散係数及び境界条件をそれぞれコンピュータに入力する第９ステップと、単結晶１４のメッシュの座標及び温度と格子間原子及び空孔の拡散係数及び境界条件に基づいて拡散方程式を解くことにより単結晶１４の冷却後の格子間原子濃度及び空孔濃度の分布を求める第１０ステップと、単結晶１４内の格子間原子濃度及び空孔濃度が等量となる分布線ＣＬ_Xを求める第１１ステップと、分布線ＣＬ_Xの変曲点の座標の最大値と最小値との差ΔＺを求める第１２ステップと、単結晶製造条件のパラメータをＰ₂からＰ_Nに順次変え第２ステップから第１２ステップを繰返して分布線ＣＬ_Xの変曲点の座標の最大値と最小値との差ΔＺを求めた後にこの差ΔＺが最も小さくなる単結晶製造条件を求める第１３ステップとを含むコンピュータを用いて単結晶の点欠陥分布のシミュレーションを行う方法である。
【０００８】
この請求項１に記載された単結晶の点欠陥分布のシミュレーション方法では、パラメータＰ₁の単結晶製造条件で融液１２の対流を考慮して単結晶１４の内部温度分布を計算により求め、この単結晶１４の内部温度分布に基づきかつ単結晶１４内の点欠陥の拡散を考慮して単結晶１４内の点欠陥分布を計算により求めた後に、単結晶１４内の格子間原子濃度及び空孔濃度が等量となる分布線ＣＬ_Xを計算により求めて、この分布線ＣＬ_Xの変曲点Ｑの最大値及び最小値の差ΔＺを計算により求める。次に上記単結晶製造条件のパラメータをＰ₂から順にＰ_Nまで変え、上記と同様に単結晶１４内の格子間原子濃度及び空孔濃度が等量となる分布線ＣＬ_Xを計算により求めて、これらの分布線ＣＬ_Xの変曲点Ｑの最大値及び最小値の差ΔＺを計算によりそれぞれ求める。更に上記分布線ＣＬ_Xの変曲点Ｑの最大値及び最小値の差ΔＺが最も小さくなる単結晶製造条件を求める。この単結晶製造条件で単結晶１４を引上げると、無欠陥領域が半径方向に均一に分布した単結晶１４を引上げることができる。
【０００９】
【発明の実施の形態】
次に本発明の実施の形態を図面に基づいて説明する。
図４に示すように、シリコン単結晶引上げ機１１のチャンバ内には、シリコン融液１２を貯留する石英るつぼ１３が設けられる。この石英るつぼ１３は図示しないが黒鉛サセプタ及び支軸を介してるつぼ駆動手段に接続され、るつぼ駆動手段は石英るつぼ１３を回転させるとともに昇降させるように構成される。また石英るつぼ１３の外周面は石英るつぼ１３から所定の間隔をあけてヒータ（図示せず）により包囲され、このヒータは保温筒（図示せず）により包囲される。ヒータは石英るつぼ１３に投入された高純度のシリコン多結晶体を加熱・溶融してシリコン融液１２にする。またチャンバの上端には図示しないが円筒状のケーシングが接続され、このケーシングには引上げ手段が設けられる。引上げ手段はシリコン単結晶１４を回転させながら引上げるように構成される。更にチャンバ内にはアルゴンガス等の不活性ガスが流通し、引上げ中のシリコン単結晶１４は不活性ガスの整流やヒータからの輻射熱の遮断等を目的とするヒートキャップ（図示せず）により包囲される。
【００１０】
このように構成されたシリコン単結晶引上げ機１１におけるシリコン単結晶１４の点欠陥分布のシミュレーション方法を図１〜図５に基づいて説明する。
先ず第１ステップとして、シリコン単結晶引上げ機１１によりシリコン単結晶１４を引上げるときの単結晶製造条件を任意に定義する。この単結晶製造条件は引上げ機１１によりシリコン単結晶１１を引上げるときに、後述する引上げ機１１のホットゾーンにフィードバックされる変数を一定の間隔で変量させたパラメータ群Ｐ₁，Ｐ₂，…，Ｐ_Nである。また単結晶製造条件としては、シリコン単結晶１４の引上げ速度、シリコン単結晶１４の回転速度、石英るつぼ１３の回転速度、アルゴンガスの流量、ヒートキャップを構成する部材の形状及び材質、ヒートキャップの下端及びシリコン融液１２表面間のギャップ、ヒータ出力等が挙げられる。
【００１１】
第２ステップとして、シリコン単結晶１４を所定長さＬ₁（例えば１００ｍｍ）まで引上げた状態における引上げ機１１のホットゾーンの各部材、即ちチャンバ，石英るつぼ１３，シリコン融液１２，シリコン単結晶１４，黒鉛サセプタ，保温筒等をメッシュ分割してモデル化する。具体的には上記ホットゾーンの各部材のメッシュ点の座標データをコンピュータに入力する。このときシリコン融液１２のメッシュのうちシリコン単結晶１４の径方向のメッシュであってかつシリコン融液１２のシリコン単結晶１４直下の一部又は全部のメッシュ（以下、径方向メッシュという。）を０．０１〜５．００ｍｍ、好ましくは０．２５〜１．００ｍｍに設定する。またシリコン融液１２のメッシュのうちシリコン単結晶１４の長手方向のメッシュであってかつシリコン融液１２の一部又は全部のメッシュ（以下、長手方向メッシュという。）を０．０１〜５．００ｍｍ、好ましくは０．１〜０．５ｍｍに設定する。
【００１２】
径方向メッシュを０．０１〜５．００ｍｍの範囲に限定したのは、０．０１ｍｍ未満では計算時間が極めて長くなり、５．００ｍｍを越えると計算が不安定になり、繰返し計算を行っても固液界面形状が一定に定まらなくなるからである。また長手方向メッシュを０．０１〜５．００ｍｍの範囲に限定したのは、０．０１ｍｍ未満では計算時間が極めて長くなり、５．００ｍｍを越えると固液界面形状の計算値が実測値と一致しなくなるからである。なお、径方向メッシュの一部を０．０１〜５．００の範囲に限定する場合には、シリコン単結晶１４直下のシリコン融液１２のうちシリコン単結晶１４外周縁近傍のシリコン融液１２を上記範囲に限定することが好ましく、長手方向メッシュの一部を０．０１〜５．００の範囲に限定する場合には、シリコン融液１２の液面近傍及び底近傍を上記範囲に限定することが好ましい。
【００１３】
第３ステップとして上記ホットゾーンの各部材毎にメッシュをまとめ、かつこのまとめられたメッシュに対して各部材の物性値をそれぞれコンピュータに入力する。例えば、チャンバがステンレス鋼にて形成されていれば、そのステンレス鋼の熱伝導率，輻射率，粘性率，体積膨張係数，密度及び比熱がコンピュータに入力される。またシリコン単結晶１４の引上げ長と、この引上げ長に対応するシリコン単結晶１４の引上げ速度等の単結晶製造条件としてパラメータＰ₁と、後述する乱流モデル式（２）の乱数パラメータＣとをコンピュータに入力する。
【００１４】
第４ステップとして、ホットゾーンの各部材の表面温度分布をヒータの発熱量及び各部材の輻射率に基づいてコンピュータを用いて計算により求める。即ち、ヒータの発熱量を任意に設定してコンピュータに入力するとともに、各部材の輻射率から各部材の表面温度分布をコンピュータを用いて計算により求める。次に第５ステップとしてホットゾーンの各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式（１）をコンピュータを用いて解くことにより各部材の内部温度分布を計算により求める。ここでは、記述を簡単にするためｘｙｚ直交座標系を用いたが、実際の計算では円筒座標系を用いる。
【００１５】
【数１】

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

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

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

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

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

式（９）及び（１０）のそれぞれ右側の第１項のＤ_i及びＤ_vは、次の式（１１）及び（１２）に示すように拡散係数を有するFickian拡散を表す。
【００２５】
【数７】

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

は熱拡散による格子間シリコン及び空孔の活性化エネルギーであり、式（９）及び式（１０）のそれぞれ右側の第３項のｋ_ivは格子間シリコン及び空孔ペアの再結合定数である。
【００２６】
第１１ステップとして、第１０ステップで求めたシリコン単結晶１４内の格子間シリコン濃度及び空孔濃度からこれらの濃度が等量となる分布線ＣＬ_X、即ちシリコン単結晶１４が冷却した後のシリコン単結晶１４内の格子間シリコン濃度及び空孔濃度が等量となる分布線ＣＬ_Xを計算により求める（図４及び図５）。その後、第１２ステップとして、上記分布線ＣＬ_Xの変曲点を計算により求める。例えば変曲点が５つある場合には（図５）、各変曲点Ｑ(1,X)、Ｑ(2,X)、Ｑ(3,X)、Ｑ(4,X)及びＱ(5,X)の座標を(ｒ_Q(1,X)，ｚ_Q(1,X))、(ｒ_Q(2,X)，ｚ_Q(2,X))、(ｒ_Q(3,X)，ｚ_Q(3,X))、(ｒ_Q(4,X)，ｚ_Q(4,X))及び(ｒ_Q(5,X)，ｚ_Q(5,X))とすると、これらの座標のうちｚ座標ｚ_Q(1,X)、ｚ_Q(2,X)、ｚ_Q(3,X)、ｚ_Q(4,X)及びｚ_Q(5,X)の最大値ｚ_Q(3,X)と最小値ｚ_Q(2,X)との差ΔＺ(X)を求める。
【００２７】
第１３ステップとして、単結晶製造条件のパラメータをＰ₂に変え第２ステップから第１２ステップを実行し、単結晶製造条件のパラメータをＰ₃に変え第２ステップから第１２ステップを実行するというように、単結晶製造条件のパラメータがＰ_Nになるまで上記第２ステップから第１２ステップを繰返した後に、分布線ＣＬ_Xの変曲点Ｑの座標の最大値と最小値との差ΔＺ(X)をそれぞれ計算により求める。更に上記各差ΔＺ(X)のうち最も小さくなる単結晶製造条件を計算により求める。
【００２８】
このように計算して求められたシリコン単結晶引上げ機１１における単結晶製造条件は、分布線ＣＬ_Xが略水平にフラットになるため、この単結晶製造条件でシリコン単結晶を引上げると、格子間原子優勢無欠陥領域［Ｐ_I］及び空孔優勢無欠陥領域［Ｐ_V］を含む無欠陥領域［Ｐ］が半径方向に均一に分布したシリコン単結晶１４を引上げることができる。換言すれば、上述のようにコンピュータを用いた計算によりシリコン単結晶１４の最適な単結晶製造条件を設定することができるので、実際にシリコン単結晶を引上げてライフタイム等を測定する実験は上記単結晶製造条件が最適なものであるか否かの確認のためにだけ行えば済む。この結果、上記実験の回数を低減できるので、比較的少ない時間と労力で、最適な単結晶製造条件を得ることができる。
なお、この実施の形態では、シリコン単結晶を挙げたが、ＧａＡｓ単結晶，ＩｎＰ単結晶，ＺｎＳ単結晶若しくはＺｎＳｅ単結晶でもよい。
【００２９】
【発明の効果】
以上述べたように、本発明によれば、所定の単結晶製造条件で融液の対流を考慮することにより、計算した結晶及び融液の固液界面形状が実際の単結晶の引上げ時の形状とほぼ一致するようにして単結晶の内部温度を求め、この単結晶の内部温度分布に基づきかつ単結晶内の点欠陥の拡散を考慮して単結晶内の点欠陥分布を求めた後に、単結晶内の格子間原子濃度及び空孔濃度が等量となる分布線を求めてこの分布線の変曲点の最大値及び最小値の差を求め、上記単結晶製造条件を変え上記と同様にして単結晶内の格子間原子濃度及び空孔濃度が等量となる分布線を求め、この分布線の変曲点の最大値及び最小値の差を求め、更に分布線の変曲点の最大値及び最小値の差が最も小さくなる単結晶製造条件を求めたので、この単結晶製造条件で単結晶を引上げれば、無欠陥領域が半径方向に均一に分布した単結晶を引上げることができる。このため単結晶の最適な引上げ条件を設定するのに数多くの実験を繰返す必要のあった従来の無欠陥の単結晶インゴットと比較して、本発明では、単結晶の最適な引上げ条件を設定するための実験の回数を低減できる。この結果、本発明では、比較的少ない時間と労力で、最適な単結晶製造条件を得ることができる。
【図面の簡単な説明】
【図１】本発明実施形態シリコン単結晶の無欠陥領域シミュレーション方法の前段を示すフローチャート。
【図２】そのシリコン単結晶の無欠陥領域シミュレーション方法の中段を示すフローチャート。
【図３】そのシリコン単結晶の無欠陥領域シミュレーション方法の後段を示すフローチャート。
【図４】本発明実施形態のシリコン融液をメッシュ構造とした引上げ機の要部断面図。
【図５】本発明実施形態のシリコン単結晶内の格子間シリコン及び空孔の分布を示す説明図。
【符号の説明】
１１シリコン単結晶引上げ機
１２シリコン融液
１４シリコン単結晶
Ｓシリコンの三重点
ＣＬ_X 分布線[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a method for computer simulation of a defect-free region of a single crystal such as silicon that is pulled by the Czochralski (hereinafter referred to as CZ) method.
[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 in Japanese Patent Laid-Open No. 11-1393 corresponding to US Pat. No. 6,045,610 Has been. 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]
[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 do not occur, the experiment of measuring the lifetime of the single crystal ingot by actually pulling it up for each type of puller and the thickness of the single crystal ingot must be repeated a number of times. 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. A one-dimensional model that simplifies Boronkov's theory of controlling V / G (mm ² / min · ° C) when the temperature gradient near the liquid interface is G (° C / mm) (only in the pulling direction of the single crystal) Therefore, the point defect diffusion (outward diffusion) in the radial direction of the single crystal is not taken into account, and the concentration gradient of the point defect in the radial direction of the single crystal cannot be obtained. It was.
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.
The object of the present invention is to reduce the number of experiments of actually pulling up a single crystal, thereby enabling the setting of optimum pulling conditions for the single crystal with relatively little time and effort, and no defects in the single crystal. It is to provide a region simulation method.
[0007]
[Means for Solving the Problems]
As shown in FIGS. 1 to 5, the invention according to claim 1 is a parameter group P in which variables fed back as a single crystal production condition when the single crystal 14 is pulled up by the puller 11 are varied at regular intervals. ₁ , P ₂ ,..., P _N are arbitrarily defined, and the hot zone of the puller 11 of the single crystal 14 in a state where the single crystal 14 is pulled up to a predetermined length by the puller 11 is modeled in a mesh structure. a second step of the parameter group P ₁ with physical properties of each member to the mesh summarized and the summary of the mesh for each member of the hot zone, P _2, ..., the parameters P ₁ of the P _N to the computer A third step of inputting, a fourth step of obtaining a surface temperature distribution of each member based on a heat generation amount of the heater and a radiation rate of each member, a surface temperature distribution of each member, and After obtaining the internal temperature distribution of each member by solving the heat conduction equation based on the conductivity, the turbulent model equation obtained by assuming that the melt 12 is turbulent and the Navier-Stokes equation are connected. The fifth step for further determining the internal temperature distribution of the melt 12 in consideration of the convection by solving, and the solid-liquid interface shape of the single crystal 14 and the melt 12 is determined according to the isotherm including the triple point S of the single crystal. The sixth step and the fourth to sixth steps 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 obtain the mesh coordinates and temperature of the single crystal 14. a seventh step of inputting data to the computer, respectively, the temperature distribution of the pull-up length of the single crystal 14 from the second step is changed stepwise in the pulling machine 11 repeatedly until the seventh step calculated monocrystal The eighth step of obtaining the coordinates and temperature of the mesh 4 and inputting these data to the computer, respectively, the mesh coordinates and temperature data of the single crystal 14, the diffusion coefficients of interstitial atoms and vacancies in the single crystal 14, and After cooling the single crystal 14, the ninth step of inputting boundary conditions to the computer, respectively, and solving the diffusion equation based on the mesh coordinates and temperature of the single crystal 14, diffusion coefficients of interstitial atoms and vacancies, and boundary conditions A tenth step for obtaining the distribution of the interstitial atom concentration and the vacancy concentration, an eleventh step for obtaining a distribution line CL _{X in} which the interstitial atom concentration and the vacancy concentration in the single crystal 14 are equal, and the distribution line CL. a twelfth step of obtaining a difference ΔZ between the maximum value and the minimum value of the coordinates of the inflection point of the _X, the parameters of the single crystal manufacturing condition from sequentially changing the second step from P ₂ to P _N the first Repeat steps using a computer including a thirteenth step for the difference ΔZ after obtaining a difference ΔZ between the maximum value and the minimum value of the coordinates of the inflection point of the distribution line CL _X Find the smallest single crystal manufacturing condition This is a method of simulating the point defect distribution of a single crystal.
[0008]
In the method for simulating the point defect distribution of a single crystal described in claim 1, the internal temperature distribution of the single crystal 14 is obtained by calculation in consideration of the convection of the melt 12 under the single crystal production conditions of the parameter P _1. After calculating the point defect distribution in the single crystal 14 based on the internal temperature distribution of the single crystal 14 and considering the diffusion of point defects in the single crystal 14, the interstitial atom concentration and vacancies in the single crystal 14 are calculated. seeking by calculation a distribution line CL _X concentration becomes equal amount, determined by calculating a difference ΔZ between the maximum value and the minimum value of the inflection point Q of the distribution line CL _X. Next, the parameters of the single crystal production conditions are changed from P ₂ to P _{N in} order, and the distribution line CL _X where the interstitial atom concentration and the vacancy concentration in the single crystal 14 are equal is obtained by calculation in the same manner as described above. The difference ΔZ between the maximum value and the minimum value of the inflection point Q of these distribution lines CL _X is obtained by calculation. Further, a single crystal production condition is obtained in which the difference ΔZ between the maximum value and the minimum value of the inflection point Q of the distribution line CL _X is minimized. When the single crystal 14 is pulled up under the single crystal manufacturing conditions, the single crystal 14 in which defect-free regions are uniformly distributed in the radial direction can be pulled up.
[0009]
DETAILED DESCRIPTION OF THE INVENTION
Next, embodiments of the present invention will be described with reference to the drawings.
As shown in FIG. 4, 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. Further, an inert gas such as argon gas circulates in the chamber, and the silicon single crystal 14 being pulled is surrounded by a heat cap (not shown) for the purpose of rectifying the inert gas and blocking radiant heat from the heater. Is done.
[0010]
A method for simulating the point defect distribution of 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, 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 ₁ , P ₂ ,... In which variables fed back to a hot zone of the pulling machine 11 to be described later are changed at regular intervals when the silicon single crystal 11 is pulled by the pulling machine 11. , P _N. 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.
[0011]
As a second step, each member of the hot zone of the puller 11 in a state where the silicon single crystal 14 is pulled up to a predetermined length L ₁ (for example, 100 mm), that is, the chamber, the quartz crucible 13, the silicon melt 12, and the silicon single crystal 14 , Graphite susceptor, insulation cylinder, etc. are modeled by dividing the mesh. 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.
[0012]
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.
[0013]
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, as a single crystal manufacturing condition such as a pulling length of the silicon single crystal 14 and a pulling speed of the silicon single crystal 14 corresponding to the pulling length, a parameter P ₁ and a random parameter C of a turbulent flow model equation (2) described later are used. Input to the computer.
[0014]
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, the xyz orthogonal coordinate system is used, but in the actual calculation, a cylindrical coordinate system is used.
[0015]
[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 , λ _y and λ _z Is 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 Navier-Stokes equations (3) to (5) obtained by assuming that there is, the internal flow velocity distribution of the silicon melt 12 is obtained by calculation using a computer.
[0016]
[Expression 2]

Here, κ _t is the turbulent thermal conductivity of the silicon melt 12, c is the specific heat of the silicon melt 12, Pr _t is the Prandtl number, ρ is the density of the silicon melt 12, and C Is a turbulent flow parameter, d is a distance from the wall of the quartz crucible 13 storing the silicon melt 12, and k is a square sum of fluctuation components with respect to the average flow velocity of the silicon melt 12.
[0017]
[Equation 3]

[0018]
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 ν _l is a molecular kinematic viscosity coefficient (physical property value) of the silicon melt 12. ν _t is a kinematic viscosity coefficient due to the effect of turbulent flow of the silicon melt 12, and F _x , F _y, and F _z are x, y, and z direction components of the body force acting on the silicon melt 12.
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 calculated using a computer. Ask for.
[0019]
[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, κ _l is the molecular thermal conductivity (physical property value), and κ _t is the turbulent heat calculated using equation (2). Conductivity.
[0020]
Next, as a sixth step, a silicon-solid interface 14 and a silicon-melt interface 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.
[0021]
Next, as an eighth step, after adding δ (for example, 50 mm) to the pulling length L ₁ of the silicon single crystal 14 and repeating the second to seventh steps, the temperature distribution in the pulling machine 11 is calculated. The coordinates and temperature of the mesh of the silicon single crystal 14 are obtained, and these data are stored in a computer. This eighth step is performed until the pulling length L ₁ of the silicon single crystal 14 reaches the length L ₂ . When the pulling length L ₁ of the silicon single crystal 14 reaches the length L ₂ , the process proceeds to the ninth step, and the coordinates and temperature data of the mesh of the silicon single crystal 14 are converted into the interstitial silicon and silicon in the silicon single crystal 14. Each is input to the computer along with the diffusion coefficient of the pores and the boundary conditions. Further, the process proceeds to the tenth step, and the concentration distribution of interstitial silicon and vacancies after cooling of the silicon single crystal 14 is calculated by solving the diffusion equation based on the diffusion coefficient and boundary conditions of the interstitial silicon and vacancies. Ask for.
[0022]
Specifically, the formula for calculating the interstitial silicon concentration C _i is given by the following formula (7), and the formula for calculating the vacancy density C _v is given by the following formula (8). In equations (7) and (8), the thermal equilibrium between interstitial silicon and vacancies is maintained at the crystal side, top and solid-liquid interface to calculate the evolution of concentration C _i and concentration C _v over time. Assume that
[0023]
[Equation 5]

Here, K ₁ and K ₂ are constants, E _i and E _v are the formation energy of interstitial silicon and vacancies, the shoulder letter e is the equilibrium amount, k is the Boltzmann constant, and T is the absolute temperature.
The above equilibrium equation is differentiated with respect to time, and becomes the following equations (9) and (10) for interstitial silicon and vacancies, respectively.
[0024]
[Formula 6]

The first terms D _i and D _v on the right side of the equations (9) and (10) respectively represent Fickian diffusion having diffusion coefficients as shown in the following equations (11) and (12).
[0025]
[Expression 7]

Here, ΔE _i and ΔE _v are activation energies of interstitial silicon and vacancies, respectively, and d _i and d _v are constants, respectively. In addition, in the second term on the right side of each of the equations (9) and (10),

Is the activation energy of interstitial silicon and vacancies due to thermal diffusion, and k _{iv in} the third term on the right side of equations (9) and (10) is the recombination constant of interstitial silicon and vacancy pairs, respectively. .
[0026]
As an eleventh step, the distribution line CL _{X in} which these concentrations are equal from the interstitial silicon concentration and the vacancy concentration in the silicon single crystal 14 obtained in the tenth step, that is, the silicon after the silicon single crystal 14 is cooled. A distribution line CL _{X in} which the interstitial silicon concentration and the vacancy concentration in the single crystal 14 are equal is obtained by calculation (FIGS. 4 and 5). Thereafter, as a twelfth step, an inflection point of the distribution line CL _X is obtained by calculation. For example, when there are five inflection points (FIG. 5), each inflection point Q (1, X), Q (2, X), Q (3, X), Q (4, X) and Q ( 5, X) coordinates (r _Q (1, X), z _Q (1, X)), (r _Q (2, X), z _Q (2, X)), (r _Q (3, X ), Z _Q (3, X)), (r _Q (4, X), z _Q (4, X)) and (r _Q (5, X), z _Q (5, X)) z-coordinate z _Q of the coordinates _{(1, X), z Q} (2, X), z Q (3, X), z Q (4, X) and z maximum z _Q of _Q (5, X) A difference ΔZ (X) between (3, X) and the minimum value z _Q (2, X) is obtained.
[0027]
As the thirteenth step, the parameter of the single crystal production condition is changed to P ₂ , the second step to the twelfth step are executed, the parameter of the single crystal production condition is changed to P _3, and the second step to the twelfth step are executed. In addition, after repeating the second step to the twelfth step until the parameter of the single crystal production condition becomes P _N , the difference ΔZ (X between the maximum value and the minimum value of the inflection point Q of the distribution line CL _X is obtained. ) By calculation. Further, the smallest single crystal production condition among the above differences ΔZ (X) is obtained by calculation.
[0028]
The single crystal manufacturing conditions in the silicon single crystal pulling machine 11 calculated in this way are such that the distribution line CL _X is flat in a substantially horizontal manner. It is possible to pull up the silicon single crystal 14 in which the defect-free region [P] including the interatomic dominant defect-free region [P _I ] and the hole dominant defect-free region [P _V ] is uniformly distributed in the radial direction. In other words, since the optimum single crystal manufacturing conditions of the silicon single crystal 14 can be set by calculation using a computer as described above, experiments for actually measuring the lifetime and the like by pulling up the silicon single crystal are as described above. It is only necessary to confirm whether or not the single crystal production conditions are optimum. As a result, since the number of experiments can be reduced, the optimum single crystal production conditions can be obtained with relatively little time and effort.
In this embodiment, a silicon single crystal is used, but a GaAs single crystal, an InP single crystal, a ZnS single crystal, or a ZnSe single crystal may be used.
[0029]
【The invention's effect】
As described above, according to the present invention, by considering the convection of the melt under predetermined single crystal production conditions, the calculated solid-liquid interface shape of the crystal and the melt is the shape when the actual single crystal is pulled up. After obtaining the point temperature distribution in the single crystal based on the internal temperature distribution of the single crystal and taking into account the diffusion of point defects in the single crystal, Obtain a distribution line in which the interstitial atom concentration and vacancy concentration in the crystal are equal, find the difference between the maximum value and the minimum value of the inflection point of this distribution line, change the above-mentioned single crystal production conditions, and perform the same as above. To obtain a distribution line in which the interstitial atom concentration and the vacancy concentration in the single crystal are equal, find the difference between the maximum value and the minimum value of the inflection point of the distribution line, and further determine the maximum of the inflection point of the distribution line. The single crystal manufacturing conditions that minimize the difference between the minimum value and the minimum value were determined. If Re pulled crystal can be pulled single crystal defect-free region is uniformly distributed in the radial direction. For this reason, in the present invention, the optimum pulling condition for the single crystal is set as compared with the conventional defect-free single crystal ingot that required many experiments to set the optimum pulling condition for the single crystal. The number of experiments can be reduced. As a result, in the present invention, it is possible to obtain optimum single crystal production conditions with relatively little time and effort.
[Brief description of the drawings]
FIG. 1 is a flowchart showing a first stage of a defect-free region simulation method for a silicon single crystal according to an embodiment of the present invention.
FIG. 2 is a flowchart showing a middle stage of the defect-free region simulation method for the silicon single crystal.
FIG. 3 is a flowchart showing the latter stage of the silicon single crystal defect-free region simulation method.
FIG. 4 is a cross-sectional view of a main part of a pulling machine having a mesh structure of silicon melt according to an embodiment of the present invention.
FIG. 5 is an explanatory view showing the distribution of interstitial silicon and vacancies in the silicon single crystal according to the embodiment of the present invention.
[Explanation of symbols]
11 Silicon single crystal pulling machine 12 Silicon melt 14 Silicon single crystal S Silicon triple point CL _X distribution line

Claims

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;
A second step of modeling a hot zone of the puller (11) of the single crystal (14) in a state where the single crystal (14) is pulled up to a predetermined length by the puller (11) with 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;
By changing the pulling length of the single crystal (14) stepwise, the temperature distribution in the pulling machine (11) is calculated repeatedly from the second step to the seventh step, and the mesh of the single crystal (14) is calculated. An eighth step of obtaining coordinates and temperature and inputting these data to the computer, respectively;
Ninth step of inputting the coordinates and temperature data of the mesh of the single crystal (14) and the diffusion coefficients and boundary conditions of interstitial atoms and vacancies and boundary conditions in the single crystal (14), respectively, to the computer;
The interstitial atoms after cooling of the single crystal (14) by solving the diffusion equation based on the coordinates and temperature of the mesh of the single crystal (14) and the diffusion coefficients and boundary conditions of the interstitial atoms and the vacancies A tenth step for determining a concentration and a distribution of the pore concentration;
An eleventh step of obtaining a distribution line (CL _X ) in which the interstitial atom concentration and the vacancy concentration in the single crystal (14) are equal;
A twelfth step of obtaining a difference (ΔZ) between the maximum value and the minimum value of the inflection point of the distribution line (CL _X );
The parameter of the single crystal production condition is sequentially changed from P ₂ to _PN, and the second step to the twelfth step are repeated, and the difference between the maximum value and the minimum value of the inflection point of the distribution line (CL _X ) is determined. A method of simulating a point defect distribution of a single crystal using a computer including a thirteenth step of determining a single crystal manufacturing condition that minimizes the difference (ΔZ) after determining (ΔZ).