(57)【特許請求の範囲】
【請求項１】
引上げ機(11)により単結晶(14)を所定長さまで引上げた状態における前記単結晶(14)の引上げ機(11)のホットゾーンをメッシュ構造でモデル化する第１ステップと、
前記ホットゾーンの各部材毎にメッシュをまとめかつこのまとめられたメッシュに対する前記各部材の物性値とともに前記単結晶(14)の引上げ長及びこの引上げ長に対応する前記単結晶(14)の引上げ速度をそれぞれコンピュータに入力する第２ステップと、
前記各部材の表面温度分布をヒータの発熱量及び前記各部材の輻射率に基づいて求める第３ステップと、
前記各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式を解くことにより前記各部材の内部温度分布を求めた後に融液(12)が乱流であると仮定して得られた乱流モデル式及びナビエ・ストークスの方程式を連結して解くことにより対流を考慮した前記融液(12)の内部温度分布を更に求める第４ステップと、
前記単結晶(14)及び前記融液(12)の固液界面形状を前記単結晶の三重点(S)を含む等温線に合せて求める第５ステップと、
前記第３ステップから前記第５ステップを前記三重点(S)が前記単結晶(14)の融点になるまで繰返し前記引上げ機(11)内の温度分布を計算して前記単結晶(14)のメッシュの座標及び温度を求めこれらのデータをそれぞれ前記コンピュータに入力する第６ステップと、
前記単結晶(14)の引上げ長を段階的に変えて前記第１ステップから前記第６ステップまでを繰返し前記引上げ機(11)内の温度分布を計算して前記単結晶(14)のメッシュの座標及び温度を求めこれらのデータをそれぞれ前記コンピュータに入力する第７ステップと、
前記単結晶(14)のメッシュの座標及び温度のデータと前記単結晶(14)内の空孔及び格子間原子の拡散係数及び境界条件をそれぞれ前記コンピュータに入力する第８ステップと、
前記単結晶(14)のメッシュの座標及び温度と前記空孔及び前記格子間原子の拡散係数及び境界条件に基づいて拡散方程式を解くことにより前記単結晶(14)の冷却後の前記空孔及び前記格子間原子の濃度分布を求める第９ステップと
を含むコンピュータを用いて単結晶の点欠陥分布のシミュレーションを行う方法であって、
前記乱流モデル式が次の式（２）で表されるｋｌ−モデル式であり、このモデル式の乱流パラメータＣとして０．４〜０．６の範囲内の任意の値が用いられた
ことを特徴とする単結晶及び融液の固液界面形状のシミュレーション方法。
【数２】
0004096499.tif 000011
ここで、κtは融液の乱流熱伝導率であり、ｃは融液の比熱であり、Ｐｒtはプラントル数であり、ρは融液の密度であり、ｄは融液を貯留するるつぼ壁からの距離であり、ｋは融液の平均流速に対する変動成分の二乗和である。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、チョクラルスキー（以下、ＣＺという。）法にて引上げられるシリコン等の単結晶の点欠陥分布をコンピュータシミュレーションする方法に関するものである。
【０００２】
【従来の技術】
従来、この種のシミュレーション方法として、図４に示すように、総合伝熱シミュレータを用いてＣＺ法によるシリコン単結晶４引上げ時の引上げ機１内のホットゾーン構造及びそのシリコン単結晶４の引上げ速度に基づいて、シリコン融液２の熱伝導率を操作することによりシリコン融液２の内部温度分布を予測し、この内部温度分布からシリコン単結晶４のメッシュの座標及び温度をそれぞれ求め、更にシリコン単結晶４内の格子間シリコン及び空孔の拡散係数及び境界条件に基づいて拡散方程式を解くことにより、上記格子間シリコン及び空孔の濃度分布をコンピュータを用いて求める方法が知られている。このシミュレーション方法では、ホットゾーンの各部材がメッシュ分割されてモデル化される。特にシリコン融液２のメッシュは計算時間を短くするために１０ｍｍ程度と比較的粗く設定される。
【０００３】
【発明が解決しようとする課題】
しかし、上記従来の格子間シリコン及び空孔の濃度分布のシミュレーション方法では、実際の引上げ機においては発生するシリコン融液の対流を考慮しておらず、またシリコン融液のメッシュが比較的粗いため、格子間シリコン及び空孔の濃度分布（図５（ｂ））が実測値（図５（ｅ））と大幅に相違する不具合があった。
本発明の目的は、単結晶内の点欠陥分布の計算値が実測値と極めて良く一致する、単結晶の点欠陥分布のシミュレーション方法を提供することにある。
【０００４】
【課題を解決するための手段】
請求項１に係る発明は、図１〜図３に示すように、引上げ機１１により単結晶１４を所定長さまで引上げた状態における単結晶１４の引上げ機１１のホットゾーンをメッシュ構造でモデル化する第１ステップと、ホットゾーンの各部材毎にメッシュをまとめかつこのまとめられたメッシュに対する各部材の物性値と単結晶１４の引上げ長及びこの引上げ長に対応する単結晶１４の引上げ速度をそれぞれコンピュータに入力する第２ステップと、各部材の表面温度分布をヒータの発熱量及び各部材の輻射率に基づいて求める第３ステップと、各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式を解くことにより各部材の内部温度分布を求めた後に融液１２が乱流であると仮定して得られた乱流モデル式及びナビエ・ストークスの方程式を連結して解くことにより対流を考慮した融液１２の内部温度分布を更に求める第４ステップと、単結晶１４及び融液１２の固液界面形状を単結晶の三重点Ｓを含む等温線に合せて求める第５ステップと、第３ステップから第５ステップを三重点Ｓが単結晶１４の融点になるまで繰返し引上げ機１１内の温度分布を計算して単結晶１４のメッシュの座標及び温度を求めこれらのデータをそれぞれコンピュータに入力する第６ステップと、単結晶１４の引上げ長を段階的に変えて第１ステップから第６ステップまでを繰返し引上げ機１１内の温度分布を計算して単結晶１４のメッシュの座標及び温度を求めこれらのデータをそれぞれコンピュータに入力する第７ステップと、単結晶１４のメッシュの座標及び温度のデータと単結晶１４内の空孔及び格子間原子の拡散係数及び境界条件をそれぞれコンピュータに入力する第８ステップと、単結晶１４のメッシュの座標及び温度と空孔及び前記格子間原子の拡散係数及び境界条件に基づいて拡散方程式を解くことにより単結晶１４の冷却後の空孔及び格子間原子の濃度分布を求める第９ステップとを含むコンピュータを用いて単結晶の点欠陥分布のシミュレーションを行う方法であって、乱流モデル式が次の式（２）で表されるｋｌ−モデル式であり、このモデル式の乱流パラメータＣとして０．４〜０．６の範囲内の任意の値が用いられたことを特徴とする。
【数２】
0004096499.tif 000002
ここで、κtは融液の乱流熱伝導率であり、ｃは融液の比熱であり、Ｐｒtはプラントル数であり、ρは融液の密度であり、ｄは融液を貯留するるつぼ壁からの距離であり、ｋは融液の平均流速に対する変動成分の二乗和である。
【０００５】
この請求項１に記載された単結晶の点欠陥分布のシミュレーション方法では、融液１２の対流を考慮して単結晶１４の内部温度分布を求め、この単結晶１４の内部温度分布に基づきかつ単結晶１４内の点欠陥の拡散を考慮して単結晶１４内の点欠陥分布を求めたので、単結晶１４内の点欠陥分布の計算値が実測値と極めて良く一致する。
【０００６】
【発明の実施の形態】
次に本発明の実施の形態を図面に基づいて説明する。
図３に示すように、シリコン単結晶引上げ機１１のチャンバ内には、シリコン融液１２を貯留する石英るつぼ１３が設けられる。この石英るつぼ１３は図示しないが黒鉛サセプタ及び支軸を介してるつぼ駆動手段に接続され、るつぼ駆動手段は石英るつぼ１３を回転させるとともに昇降させるように構成される。また石英るつぼ１３の外周面は石英るつぼ１３から所定の間隔をあけてヒータ（図示せず）により包囲され、このヒータは保温筒（図示せず）により包囲される。ヒータは石英るつぼ１３に投入された高純度のシリコン多結晶体を加熱・溶融してシリコン融液１２にする。またチャンバの上端には図示しないが円筒状のケーシングが接続され、このケーシングには引上げ手段が設けられる。引上げ手段はシリコン単結晶１４を回転させながら引上げるように構成される。
【０００７】
このように構成されたシリコン単結晶引上げ機１１におけるシリコン単結晶１４の点欠陥分布のシミュレーション方法を図１〜図３に基づいて説明する。
先ず第１ステップとして、シリコン単結晶１４を所定長さＬ1（例えば１００ｍｍ）まで引上げた状態におけるシリコン単結晶引上げ機１１のホットゾーンの各部材、即ちチャンバ，石英るつぼ１３，シリコン融液１２，シリコン単結晶１４，黒鉛サセプタ，保温筒等をメッシュ分割してモデル化する。具体的には上記ホットゾーンの各部材のメッシュ点の座標データをコンピュータに入力する。このときシリコン融液１２のメッシュのうちシリコン単結晶１４の径方向のメッシュであってかつシリコン融液１２のシリコン単結晶１４直下の一部又は全部のメッシュ（以下、径方向メッシュという。）を０．０１〜５．００ｍｍ、好ましくは０．２５〜１．００ｍｍに設定する。またシリコン融液１２のメッシュのうちシリコン単結晶１４の長手方向のメッシュであってかつシリコン融液１２の一部又は全部のメッシュ（以下、長手方向メッシュという。）を０．０１〜５．００ｍｍ、好ましくは０．１〜０．５ｍｍに設定する。
【０００８】
径方向メッシュを０．０１〜５．００ｍｍの範囲に限定したのは、０．０１ｍｍ未満では計算時間が極めて長くなり、５．００ｍｍを越えると計算が不安定になり、繰返し計算を行っても固液界面形状が一定に定まらなくなるからである。また長手方向メッシュを０．０１〜５．００ｍｍの範囲に限定したのは、０．０１ｍｍ未満では計算時間が極めて長くなり、５．００ｍｍを越えると固液界面形状の計算値が実測値と一致しなくなるからである。なお、径方向メッシュの一部を０．０１〜５．００の範囲に限定する場合には、シリコン単結晶１４直下のシリコン融液１２のうちシリコン単結晶１４外周縁近傍のシリコン融液１２を上記範囲に限定することが好ましく、長手方向メッシュの一部を０．０１〜５．００の範囲に限定する場合には、シリコン融液１２の液面近傍及び底近傍を上記範囲に限定することが好ましい。
【０００９】
第２ステップとして上記ホットゾーンの各部材毎にメッシュをまとめ、かつこのまとめられたメッシュに対して各部材の物性値をそれぞれコンピュータに入力する。例えば、チャンバがステンレス鋼にて形成されていれば、そのステンレス鋼の熱伝導率，輻射率，粘性率，体積膨張係数，密度及び比熱がコンピュータに入力される。またシリコン単結晶１４の引上げ長及びこの引上げ長に対応するシリコン単結晶１４の引上げ速度と、後述する乱流モデル式（２）の乱流パラメータＣとをコンピュータに入力する。
【００１０】
第３ステップとして、ホットゾーンの各部材の表面温度分布をヒータの発熱量及び各部材の輻射率に基づいてコンピュータを用いて求める。即ち、ヒータの発熱量を任意に設定してコンピュータに入力するとともに、各部材の輻射率から各部材の表面温度分布をコンピュータを用いて求める。次に第４ステップとしてホットゾーンの各部材の表面温度分布及び熱伝導率に基づいて熱伝導方程式（１）をコンピュータを用いて解くことにより各部材の内部温度分布を求める。ここでは、記述を簡単にするためｘｙｚ直交座標系を用いたが、実際の計算では円筒座標系を用いる。
【００１１】
【数１】
0004096499.tif 000003
ここで、ρは各部材の密度であり、ｃは各部材の比熱であり、Ｔは各部材の各メッシュ点での絶対温度であり、ｔは時間であり、λx，λy及びλzは各部材の熱伝導率のｘ，ｙ及びｚ方向成分であり、ｑはヒータの発熱量である。
一方、シリコン融液１２に関しては、上記熱伝導方程式（１）でシリコン融液１２の内部温度分布を求めた後に、このシリコン融液１２の内部温度分布に基づき、シリコン融液１２が乱流であると仮定して得られた乱流モデル式（２）及びナビエ・ストークスの方程式（３）〜（５）を連結して、シリコン融液１２の内部流速分布をコンピュータを用いて求める。
【００１２】
【数２】
0004096499.tif 000004
ここで、κtはシリコン融液１２の乱流熱伝導率であり、ｃはシリコン融液１２の比熱であり、Ｐｒtはプラントル数であり、ρはシリコン融液１２の密度であり、Ｃは乱流パラメータであり、ｄはシリコン融液１２を貯留する石英るつぼ１３壁からの距離であり、ｋはシリコン融液１２の平均流速に対する変動成分の二乗和である。
【００１３】
【数３】
0004096499.tif 000005
【００１４】
ここで、ｕ，ｖ及びｗはシリコン融液１２の各メッシュ点での流速のｘ，ｙ及びｚ方向成分であり、νlはシリコン融液１２の分子動粘性係数（物性値）であり、νtはシリコン融液１２の乱流の効果による動粘性係数であり、Ｆx，Ｆy及びＦzはシリコン融液１２に作用する体積力のｘ，ｙ及びｚ方向成分である。
上記乱流モデル式（２）はｋｌ（ケイエル）−モデル式と呼ばれ、このモデル式の乱流パラメータＣは０．４〜０．６の範囲内の任意の値が用いられることが好ましい。乱流パラメータＣを０．４〜０．６の範囲に限定したのは、０．４未満又は０．６を越えると計算により求めた界面形状が実測値と一致しないという不具合があるからである。また上記ナビエ・ストークスの方程式（３）〜（５）はシリコン融液１２が非圧縮性であって粘度が一定である流体としたときの運動方程式である。
上記求められたシリコン融液１２の内部流速分布に基づいて熱エネルギ方程式（６）を解くことにより、シリコン融液１２の対流を考慮したシリコン融液１２の内部温度分布をコンピュータを用いて更に求める。
【００１５】
【数４】
0004096499.tif 000006
ここで、ｕ，ｖ及びｗはシリコン融液１２の各メッシュ点での流速のｘ，ｙ及びｚ方向成分であり、Ｔはシリコン融液１２の各メッシュ点での絶対温度であり、ρはシリコン融液１２の密度であり、ｃはシリコン融液１２の比熱であり、κlは分子熱伝導率（物性値）であり、κtは式（２）を用いて計算される乱流熱伝導率である。
【００１６】
次いで第５ステップとして、シリコン単結晶１４及びシリコン融液１２の固液界面形状を図２の点Ｓで示すシリコンの三重点Ｓ（固体と液体と気体の三重点(tri-junction)）を含む等温線に合せてコンピュータを用いて求める。第６ステップとして、コンピュータに入力するヒータの発熱量を変更し（次第に増大し）、上記第３ステップから第５ステップを三重点がシリコン単結晶１４の融点になるまで繰返した後に、引上げ機１１内の温度分布を計算してシリコン単結晶のメッシュの座標及び温度を求め、これらのデータをコンピュータに記憶させる。
【００１７】
次に第７ステップとして、シリコン単結晶１４の引上げ長Ｌ1にδ（例えば５０ｍｍ）だけ加えて上記第１ステップから第６ステップまでを繰返した後に、引上げ機１１内の温度分布を計算してシリコン単結晶１４のメッシュの座標及び温度を求め、これらのデータをコンピュータに記憶させる。この第７ステップはシリコン単結晶１４の引上げ長Ｌ1が長さＬ2に達するまで行われる。シリコン単結晶１４の引上げ長Ｌ1が長さＬ2に達すると、第８ステップに移行して、シリコン単結晶１４のメッシュの座標及び温度のデータを、シリコン単結晶１４内の格子間シリコン及び空孔の拡散係数及び境界条件とともにそれぞれコンピュータに入力する。更にこれらの格子間シリコン及び空孔の拡散係数及び境界条件に基づいて拡散方程式を解くことによりシリコン単結晶１４の冷却後の格子間シリコン及び空孔の濃度分布を求める。
【００１８】
具体的には、格子間シリコンの濃度Ｃiの計算式が次の式（７）で、空孔の濃度Ｃvの計算式が次の式（８）で示される。式（７）及び式（８）において、濃度Ｃi及び濃度Ｃvの経時的進展を計算するために、格子間シリコンと空孔の熱平衡が結晶の側面、上面及び固液界面では維持されると仮定する。
【００１９】
【数５】
0004096499.tif 000007
ここで、Ｋ1及びＫ2は定数、Ｅi及びＥvはそれぞれ格子間シリコン及び空孔の形成エネルギー、肩付き文字ｅは平衡量、ｋはボルツマン定数、Ｔは絶対温度を意味する。
上記平衡式は時間で微分され、格子間シリコン及び空孔に対してそれぞれ次の式（９）及び式（１０）になる。
【００２０】
【数６】
0004096499.tif 000008
式（９）及び（１０）のそれぞれ右側の第１項のＤi及びＤvは、次の式（１１）及び（１２）に示すように拡散係数を有するFickian拡散を表す。
【００２１】
【数７】
0004096499.tif 000009
ここで△Ｅi及び△Ｅvはそれぞれ格子間シリコン及び空孔の活性化エネルギーであり、ｄi及びｄvはそれぞれ定数である。また式（９）及び式（１０）のそれぞれ右側の第２項の
【００２２】
【数８】
0004096499.tif 000010
は熱拡散による格子間シリコン及び空孔の活性化エネルギーであり、式（９）及び式（１０）のそれぞれ右側の第３項のｋivは格子間シリコン及び空孔ペアの再結合定数である。
【００２３】
このように計算して得られたシリコン単結晶１４の点欠陥分布は実測値とほぼ一致する。この結果、引上げ機１１の設計段階でこの引上げ機１１にて引上げられるシリコン単結晶１４内の点欠陥分布を予測でき、逆に引上げられるシリコン単結晶１４内の点欠陥を所望の分布にするために、引上げ機１１の設計段階で構造を検討することができる。
なお、この実施の形態では、シリコン単結晶を挙げたが、ＧａＡｓ単結晶，ＩｎＰ単結晶，ＺｎＳ単結晶若しくはＺｎＳｅ単結晶でもよい。
【００２４】
【実施例】
次に本発明の実施例を比較例とともに詳しく説明する。
＜実施例１＞
図３に示すように、石英るつぼ１３に貯留されたシリコン融液１２から直径６インチのシリコン単結晶１４を引上げる場合の、シリコン単結晶１４内の点欠陥分布を、図１及び図２のフローチャートに基づくシミュレーション方法により求めた。即ち、シリコン単結晶引上げ機１１のホットゾーンをメッシュ構造でモデル化した。ここで、シリコン融液１２のシリコン単結晶１４直下のシリコン単結晶１４の径方向のメッシュを０．７５ｍｍに設定し、シリコン融液１２のシリコン単結晶１４直下以外のシリコン単結晶１４の径方向のメッシュを１〜５ｍｍに設定した。またシリコン融液１２のシリコン単結晶１４の長手方向のメッシュを０．２５〜５ｍｍに設定し、乱流モデル式の乱流パラメータＣとして０．４５を用いた。このような条件下で、シリコン融液１２の対流を考慮してシリコン単結晶１４の内部温度分布を求め、このシリコン単結晶１４の内部温度分布に基づきかつシリコン単結晶１４内の点欠陥の拡散を考慮してシリコン単結晶１４内の点欠陥分布を求めた。
【００２５】
＜比較例１＞
図４に示すように、石英るつぼ３に貯留されたシリコン融液２から直径６インチのシリコン単結晶４を引上げる場合の、シリコン単結晶４内の点欠陥分布を従来のシミュレーション方法により求めた。即ち、シリコン単結晶引上げ機１のホットゾーンをメッシュ構造でモデル化した。ここで、シリコン融液２のシリコン単結晶４の径方向のメッシュを１０ｍｍに設定し、シリコン融液２のシリコン単結晶４の長手方向のメッシュを１０ｍｍに設定した。またシリコン融液２の対流を考慮しなかった（乱流モデル式及びナビエ・ストークスの方程式を連結した式は用いなかった。）。上記以外は実施例１と同様にコンピュータを用いてシミュレーションを行った。
【００２６】
＜比較例２＞
シリコン融液の対流を考慮したけれども、シリコン単結晶内の点欠陥の拡散を考慮しなかったことを除いて、実施例１と同様にしてコンピュータを用いてシミュレーションを行った。
＜比較例３＞
シリコン融液の対流及びシリコン単結晶内の点欠陥の拡散のいずれも考慮しなかったことを除いて、実施例１と同様にしてコンピュータを用いてシミュレーションを行った。
【００２７】
＜比較試験及び評価＞
実施例１及び比較例１〜３のシミュレーション方法によりシリコン単結晶の点欠陥分布を求めた。その結果を図５（ａ）〜（ｄ）にシリコン単結晶の点欠陥分布の実測値（図５（ｅ））とともに示す。
図５から明らかなように、比較例１〜３のシミュレーション方法で得られたシリコン単結晶の点欠陥分布（図５（ｂ）〜（ｄ））は実測値（図５（ｅ））と大幅に相違しているのに対し、実施例１のシミュレーション方法で得られたシリコン単結晶の点欠陥分布（図５（ａ））は実測値とほぼ一致していることが判った。
【００２８】
【発明の効果】
以上述べたように、本発明によれば、融液の対流を考慮して単結晶の内部温度分布を求め、この単結晶の内部温度分布に基づきかつ単結晶内の点欠陥の拡散を考慮して単結晶内の点欠陥分布を求めたので、単結晶内の点欠陥分布の計算値が実測値と極めて良く一致する。この結果、単結晶引上げ機の設計段階でこの引上げ機にて引上げられる単結晶内の点欠陥分布を予測でき、逆に引上げられる単結晶内の点欠陥を所望の分布にするために、引上げ機の設計段階で構造を検討することができる。
【図面の簡単な説明】
【図１】本発明実施形態シリコン単結晶の点欠陥分布のシミュレーション方法の前半を示すフローチャート。
【図２】そのシリコン単結晶の点欠陥分布のシミュレーション方法の後半を示すフローチャート。
【図３】本発明のシリコン融液をメッシュ構造としたシリコン単結晶の引上げ機の要部断面図。
【図４】従来例のシリコン融液をメッシュ構造としたシリコン単結晶の引上げ機の要部断面図。
【図５】実施例１、比較例１（従来例）、比較例２、比較例３及び実際に測定したシリコン単結晶内の格子間シリコン及び空孔の分布を示す縦断面図。
【符号の説明】
１１  シリコン単結晶引上げ機
１２  シリコン融液
１４  シリコン単結晶
Ｓ  シリコンの三重点
【図１】
0004096499.tif 000012
【図２】
0004096499.tif 000013
【図３】
0004096499.tif 000014
【図４】
0004096499.tif 000015
【図５】
0004096499.tif 000016