"レーザースパイクアニールのレーザービームスキャンを受ける100を主面とするシリコン基板を複数のメッシュ要素に細分化してメッシュ化する第1ステップと、前記レーザービームのスキャン速度、前記レーザービームのサイズ、前記レーザービームのエネルギー及びスキャン経路から、前記レーザービームが照射されるメッシュ要素のエネルギーの時間変化を計算する第2ステップと、前記各メッシュ要素の座標データ、物性値、前記基板の初期温度及び前記第2ステップで計算されたエネルギーの時間変化をそれぞれコンピュータに入力する第3ステップと、前記基板の初期温度、熱伝導率に基づいたレーザービームの入射エネルギーを負荷として、前記基板の温度分布の時間変化を求める第4ステップと、前記第4ステップで求められた温度分布を利用して応力計算を実施し、前記基板に生じる熱応力分布の時間変化を求める第5ステップと、前記メッシュ化された基板の任意の位置Aを設定し、前記第4ステップで求められた温度分布及び前記第5ステップの基板に生じた熱応力分布から、設定した任意の位置Aにおける深さ方向の分布を各時間で選び出す第6ステップと、前記第6ステップで選び出した温度分布と前記基板中に含まれる酸素析出物サイズ及び残存酸素濃度から、酸素析出物から発生する転位の臨界応力を求める第7ステップと、前記第6ステップで選び出した熱応力分布の応力成分σxx、σyy、τxyから、シリコン単結晶の主スリップ面である111面の[110]方向のせん断応力成分を計算する第8ステップと、前記第8ステップのせん断応力と前記第7ステップの臨界応力とを各深さ位置ごとに比較し、比較結果から各深さ位置の転位長さを求める第9ステップと、前記第9ステップで求めた転位長さから、スリップ発生の深さを計算する第10ステップとを含むコンピュータを用いたシミュレーションによりレーザースパイクアニールを施す際に生じる酸素析出物から発生する転位を予測する方法。請求項1記載のシミュレーションによる予測方法で得られる予測結果に基づいて、酸素析出物から発生する転位の発生を最小限にすることが可能な条件をレーザースパイクアニール処理条件として選定することを特徴とするレーザースパイクアニール処理条件の選定方法。請求項1記載のシミュレーションによる予測方法で得られる予測結果に基づいて、酸素析出物から発生する転位が基板表面に達しない範囲の酸素析出核密度、並びに酸素濃度を有するシリコン基板をレーザースパイクアニール処理対象基板として選定することを特徴とするシリコン基板の選定方法。請求項1記載のシミュレーションによる予測方法で得られる予測結果に基づいて、酸素析出物から発生する転位の基板表面に達する密度が小さい範囲の酸素析出核密度、並びに酸素濃度を有するシリコン基板をレーザースパイクアニール処理対象基板として選定することを特徴とするシリコン基板の選定方法。このような問題に対処するためには、LSA工程中にシリコン基板でどのような物理的現象が起きているのかを把握し、そのメカニズムを理解しなければならない。LSA工程中の実際の現象は局所的かつ瞬間的なので、測定上の問題などで温度分布、応力分布などの実験的なデータを得るのはほぼ不可能に近い。このような観点から、実際に行われるLSA工程と同様にコンピュータ上でシミュレーションできるモデルを開発し、これを再現することにより、その中で起きる現象を把握するのが重要とされる。しかしながら、これまでLSA工程をシミュレーションできる技術は知られていなかった。本発明の目的は、LSA工程でのシリコン基板内に生じる物理的現象を正確に予測できる、シミュレーションによるLSAを施す際に生じるBMDからの発生転位予測方法を提供することにある。本発明の第1の観点は、図1に示すように、レーザースパイクアニールのレーザービームスキャンを受ける100を主面とするシリコン基板を複数のメッシュ要素に細分化してメッシュ化する第1ステップと、レーザービームのスキャン速度、レーザービームのサイズ、レーザービームのエネルギー及びスキャン経路から、レーザービームが照射されるメッシュ要素のエネルギーの時間変化を計算する第2ステップと、各メッシュ要素の座標データ、物性値、基板の初期温度及び上記第2ステップで計算されたエネルギーの時間変化をそれぞれコンピュータに入力する第3ステップと、基板の初期温度、熱伝導率に基づいたレーザービームの入射エネルギーを負荷として、基板の温度分布の時間変化を求める第4ステップと、上記第4ステップで求められた温度分布を利用して応力計算を実施し、基板に生じる熱応力分布の時間変化を求める第5ステップと、メッシュ化された基板の任意の位置Aを設定し、上記第4ステップで求められた温度分布及び上記第5ステップの基板に生じた熱応力分布から、設定した任意の位置Aにおける深さ方向の分布を各時間で選び出す第6ステップと、図2に示すように、上記第6ステップで選び出した温度分布と基板中に含まれるBMDサイズ及び残存酸素濃度から、BMDから発生する転位の臨界応力を求める第7ステップと、上記第6ステップで選び出した熱応力分布の応力成分σxx、σyy、τxyから、シリコン単結晶の主スリップ面である111面の[110]方向のせん断応力成分を計算する第8ステップと、上記第8ステップのせん断応力と上記第7ステップの臨界応力とを各深さ位置ごとに比較し、比較結果から各深さ位置の転位長さを求める第9ステップと、上記第9ステップで求めた転位長さから、スリップ発生の深さを計算する第10ステップとを含むコンピュータを用いたシミュレーションによりレーザースパイクアニールを施す際に生じるBMDから発生する転位を予測する方法である。本発明の第2の観点は、第1の観点に基づくシミュレーションによる予測方法で得られる予測結果に基づいて、BMDから発生する転位の発生を最小限にすることが可能な条件をレーザースパイクアニール処理条件として選定することを特徴とするレーザースパイクアニール処理条件の選定方法である。本発明の第3の観点は、第1の観点に基づくシミュレーションによる予測方法で得られる予測結果に基づいて、BMDから発生する転位が基板表面に達しない範囲の酸素析出核密度、並びに酸素濃度を有するシリコン基板をレーザースパイクアニール処理対象基板として選定することを特徴とするシリコン基板の選定方法である。本発明の第4の観点は、第1の観点に基づくシミュレーションによる予測方法で得られる予測結果に基づいて、BMDから発生する転位の基板表面に達する密度が小さい範囲の酸素析出核密度、並びに酸素濃度を有するシリコン基板をレーザースパイクアニール処理対象基板として選定することを特徴とするシリコン基板の選定方法である。本発明のシミュレーションによるLSAを施す際に生じる酸素析出物からの発生転位予測方法では、第1ステップから第10ステップまでの各ステップを経ることにより、LSA工程中にシリコン基板内に生じる温度分布や応力分布を正確に計算し、計算された温度分布や応力分布を用いてLSA工程中にシリコン基板内に生じる物理的現象、即ちシリコン基板に含まれるBMDから発生する転位の長さを正確に予測することができる。この予測結果に基づいて、LSA工程中に発生するシリコン基板の割れやBMDから発生する転位の発生によるシリコン基板の反りといった不具合に対処することができる。具体的には、このシミュレーションを用いて得られた予測結果に基づき、BMDから発生する転位の発生を最小限にすることが可能な条件をLSA処理条件として選定することで、LSAによる転位の発生を抑制することができる。また、このシミュレーションを用いて得られた予測結果に基づき、BMDから発生する転位が基板表面に達しない範囲のBMD核密度、並びに酸素濃度を有するシリコン基板をLSA処理対象基板として選定することで、LSAによる転位の発生を抑制することができる。更に、このシミュレーションを用いて得られた予測結果に基づき、BMDから発生する転位の基板表面に達する密度が小さい範囲の、BMD核密度、並びに酸素濃度を有するシリコン基板をLSA処理対象基板として選定することで、LSAによる転位の発生を抑制することができる。次に本発明を実施するための形態を図面に基づいて説明する。本発明は、半導体デバイス製造工程の一つであるLSA工程をコンピュータを用いたシミュレーションにより、LSA工程でのシリコン基板内に生じる温度分布や応力分布を正確に計算し、計算された温度分布や応力分布を用いてLSA工程中に生じるシリコン基板に含まれるBMDから発生する転位を予測することを特徴とする。これまでLSA工程をシミュレーションできる、レーザービームを利用したアニールの具体的な計算方法が実現しなかった。また、これまで光をシリコン基板に照射するプロセスのシミュレーションでは常温での光学定数が利用されているのが通常であった。LSA工程で使用されるCO2レーザーの波長である10.6μmでは、シリコンでの光の吸収は自由キャリア吸収FreeCarrierAbsorptionによるため、常温ではシリコン基板を透過するが、温度が上がるとそれにより熱的に発生する自由キャリア電子ホールによって吸収が増加し、光の侵入深さは温度が上がるにつれて短くなる。このため、この現象を取り入れて計算に適用しないと正確な計算ができない。本発明では、温度変化によって発生するキャリアの濃度を考慮し、レーザービームの侵入深さを計算し、レーザービームに該当するエネルギーを深さ方向に与えることで、より物理的に妥当な計算を可能とした。以下、各ステップごとに詳述する。なお、本発明のシミュレーションモデルは、有限要素法などによる汎用構造解析ソフトなど、一般に用いられる計算機シミュレーションにより行う。具体的な有限要素法などによる汎用構造解析ソフトとしては、汎用有限要素法シミュレーターであるABAQUSABAQUS社、ANSYSANSYS社、GENESISVanderplaatsResearchDevelopment社、JMAGJMAG社等が例示される。第1ステップ先ず、図1に示すように、LSAのレーザービームスキャンを受ける100を主面とするシリコン基板を複数のメッシュ要素に細分化してメッシュ化する。メッシュ化するシリコン基板の計算範囲の高さ及び幅は、シリコン基板の厚み、レーザービームサイズ、スキャン速度、シミュレーション時間、コンピュータの計算能力などを考慮して決定する。細分化するメッシュ要素は小さければ小さいほど、LSA工程で実際に起きる現象に近い正確な計算を行うことができるが、その一方で、計算するメッシュ要素の数が多くなるとコンピュータによる処理が遅くなるので、正確な結果を効率的に得られるような大きさを決定する必要がある。複数のメッシュ要素に細分化された計算範囲のX軸の方向及びY軸の方向は、任意であって、どの方向でも設定可能である。なお、図4では、一例として、計算範囲のX軸の方向を[110]方向、Y軸の方向を[100]方向と一致する面と仮定した。Y軸方向を[100]方向と一致する面と仮定したのは、100基板が一般的なシリコン基板であるためである。また、X軸方向を[110]方向と一致する面と仮定したのは、便利上のもので一般的に100基板にした場合、[110]ノッチが下になるので、横スキャンすると[110]方向スキャンになるためである。レーザービームを基板表面に照射する動作を模擬するためには、横方向に対して一定の間隔を維持するのが便利なので、各メッシュ要素の幅は一定とし、また、表面部をより詳細に計算するために、表面に近づくほどメッシュ要素の高さを小さくすることが好ましい。各メッシュ要素の幅がレーザービーム幅よりも十分に小さくなるようにメッシュ化を行うことが好ましい。また、各メッシュ要素の幅は、レーザービーム幅に対し、倍数の逆数になるように、その大きさを調整することが好ましい。例えば、図5に示すように、レーザービーム4のビームサイズが120μmである場合、メッシュ要素8個分の幅がビームサイズに相当するように、メッシュ要素1個の幅を15μmに設定する。正確な結果を得るためには、少なくともメッシュ要素5個分以上の幅がビーム幅に相当するように、メッシュ要素の幅を設定することが好ましい。第2ステップ次いで、図1に戻って、レーザービームのスキャン速度、レーザービームのサイズ、レーザービームのエネルギー及びスキャン経路から、レーザービームが照射されるメッシュ要素のエネルギーの時間変化を計算する。レーザービームを基板表面に照射する動作を模擬するためには、レーザービームの特性である小さい領域ビームサイズ領域でエネルギーを局所的に与えること、そして、レーザービームが入射されて熱に変換するのを表現する必要がある。本発明では、それをエネルギー流速の形態で表面に或いはある深さまでに与える。具体的には、領域の表面に与えるときは表面熱流速J/cm2s、領域全体に与えるときは物体熱流速J/cm3sで与える。そして、図5に示すようなレーザービームの幅に相当する複数のメッシュ要素の表面部だけに負荷として表面熱流速J/cm2sを与えれば、レーザービームがその部分に当たる効果を実現できる。物体熱流速J/cm3sで与える場合は、メッシュ要素の体積分にエネルギーを与えて熱源とすることで、レーザービームがその部分に当たる効果を実現できる。状況に応じて、表面熱流速か物体熱流速を表面部に対して与える。次に、レーザービームのスキャン動作を実現するためには、スキャン動作にあわせて、各位置でのエルギーの時間変化を与える必要がある。本発明では、図6に示すような時間変化度を利用して、各メッシュ要素への負荷の時間変化を求める。具体的には、表面熱流速をa、時間変化度をbとするとき、表面熱流速aと時間変化度bとの積、即ち、表面熱流速a×時間変化度bを求めることで、各メッシュ要素への負荷の時間変化を求める。これにより、図7に示すような、レーザービームが照射されるメッシュ要素のエネルギーの時間変化を表現できる。第3ステップ次に、図1に戻って、各メッシュ要素の座標データ、物性値、基板の初期温度及び上記第2ステップで計算されたエネルギーの時間変化をそれぞれコンピュータに入力する。物性値は、シリコンの密度2.3g/cm3、ヤング率1.7GPa、ポアソン比0.26、熱伝導度、比熱、熱膨張係数を入力する。本発明では、前述したように、温度変化によって発生するキャリアの濃度を考慮するため、熱伝導度、比熱、熱膨張係数は、図8~図10に示す温度依存性データを使用する。第4ステップ次に、図1に戻って、基板の初期温度、熱伝導率に基づいたレーザービームの入射エネルギーを負荷として、基板の温度分布の時間変化を求める。基板の温度分布の時間変化は、次の式1に示す二次元x,yに対しての熱伝導方程式を解くことにより求められ、例えば、図11に示すような温度分布の時間変化や、図12に示すようなメッシュ化された基板の温度分布が得られる。上記式1中、ρはシリコンの密度、cはシリコンの比熱、kx、kyは熱伝導度のx,y成分、Tは各メッシュ要素での温度、tは時間である。第5ステップ次に、図1に戻って、上記第4ステップで求められた温度分布を利用して応力計算を実施し、基板に生じる熱応力分布の時間変化を求める。この第5ステップでは、先ず、次の式2により各メッシュ要素での熱歪みεthを求める。そして、次の式3に示す、熱歪みεthとヤング率ETとの積から各メッシュ要素での応力σが求められる。この応力σから、X方向引張応力成分σxx、Y方向引張応力成分σyy、X方向又はY方向せん断応力成分τxyが算出される。上記関係式を解くことで、例えば、図13a~図13cに示すようなメッシュ化された基板の応力分布が得られる。上記式2中、αTは熱膨張係数、Tは現在の温度、TIは初期温度、T0は熱膨張係数に対する基準温度であり、この基準温度での熱膨張はゼロと仮定される。第6ステップ次に、図1に戻って、メッシュ化された基板の任意の位置Aを設定し、上記第4ステップで求められた温度分布及び上記第5ステップの基板に生じた熱応力分布から、設定した任意の位置Aにおける深さ方向の分布を各時間で選び出す。この第6ステップでは、例えば図14に示すように、メッシュ化された基板のある任意の位置Aを設定する。そして、この設定した任意位置Aにおける深さ方向の分布を各時間で選び出す。このように任意位置Aにおける深さ方向の分布を選び出すことで、例えば、図15に示すような温度分布などが得られる。選び出した温度分布は、臨界応力BMDからの転位を計算するために用いられる。また、選び出した応力分布は、臨界応力との比較に用いられる。第7ステップ次に、図2に示すように、上記第6ステップで選び出した温度分布と基板中に含まれるBMDサイズ及び残存酸素濃度から、BMDから発生する転位の臨界応力を求める。この第7ステップでは、上記第6ステップで選び出した温度分布と基板のBMDサイズと残存酸素濃度とをフィッティング式に入力して、BMDから発生する転位の臨界応力τBMDを求める。なお、臨界応力とは、あるBMDから転位が発生し始める応力をいう。次の式4は、本出願人が行った図16に示すような実験結果を基に、板状酸素析出物に対して求められた関連フィッティング式であり、本発明で想定される熱処理の状況から、次の式4に示す板状酸素析出物に対するフィッティング式を用いることは妥当である。この関連フィッティング式に上記第6ステップで選び出した温度分布と、基板のBMDサイズと残存酸素濃度とを入力することで、例えば、図17に示すようなBMDから発生する転位の臨界応力τBMDを求めることができる。上記式4中、Coは残存酸素濃度であり、1.1×1018cm3oldASTMと仮定する。また、LはBMDサイズであり、100nmと仮定する。kはボルツマン定数であり、Tは温度分布であり、絶対温度Kを単位とする。第8ステップ次に、図2に戻って、上記第6ステップで選び出した熱応力分布の応力成分σxx、σyy、τxyから、シリコン単結晶の主スリップ面である111面の[110]方向のせん断応力成分τxyを計算する。この第8ステップでは、二次元x,yで計算された応力成分σxx、σxy、τxyから、111面の[110]方向のせん断応力τxyを計算する。計算方法は次に示す通りである。x-y平面でのせん断応力成分はτxyだが、x方向を[110]、y方向を[100]方向としたら、この面と111面が作る角度は54.7度なので、次の式5に示す関係が導かれる。そしてこの面で110方向は、更に30°の方向になるので、次の式6に示す関係が導かれる。従って、111面の[110]方向のせん断応力はτxyとなる。よって、上記第6ステップで選び出した熱応力分布の応力成分σxx、σyy、τxyから、式5及び式6に示す関係式に基づき、せん断応力τxyを求めることができる。第9ステップ次に、図2に戻って、上記第8ステップのせん断応力τxyと上記第7ステップの臨界応力τBMDとを各深さ位置ごとに比較し、比較結果から各深さ位置の転位長さを求める。この第9ステップでは、例えば、図18に示すように、各深さ位置におけるせん断応力τxyと臨界応力τBMDとを時間ごとに比較し、比較結果から例えば、図19に示すような各深さ位置の転位長さを求める。また転位の増殖は、せん断応力τxyが臨界応力τBMDより大きくなると転位が発生し、せん断応力τxyが臨界応力τBMDより大きい限り転位の移動増殖が続き、せん断応力τxyが臨界応力τBMDより小さくなると転位は停止するものと仮定する。比較方法は具体的には次に示す手法により行われる。図3に示すように、先ず、ある深さ位置におけるせん断応力τxyと臨界応力τBMDとを時間ごとに比較し、その結果がせん断応力τxy>臨界応力τBMDの関係を満たさない場合、t0=t0Δtとして時間を進め、進めた時間におけるせん断応力τxyと臨界応力τBMDとを比較する。一方、せん断応力τxy>臨界応力τBMDの関係を満たす場合、次の式7に示す転位速度式から計算されたΔLを転位長さLdに加算していく。なお、式7に示す転位速度式から求められた転位速度VとΔtとの積から時間当たりの転位移動距離ΔLが求まる。次に、t0>スキャン時間の関係を満たさない場合、t0=t0Δtとして時間を進め、進めた時間におけるせん断応力τxyと臨界応力τBMDとを比較する。そして、t0>スキャン時間の関係を満たす場合、ある深さ位置の比較を終了する。この時点までの計算された転位長さLdをある深さ位置における転位長さとする。そして、残りの他の深さ位置についての比較も上記比較方法と同様にして行い、残りの他の深さ位置における転位長さを求める。第10ステップ次に、図2に戻って、上記第9ステップで求めた転位長さから、スリップ発生の深さを計算する。上記第9ステップで求めた各深さ位置の転位長さから、図19に示すような転位長さ曲線が得られる。この転位長さ曲線から、ある深さ位置にBMDが存在するときの、このBMDから生じる転位の長さを予測できる。そして、転位長さLd>深さDの関係から、スリップ発生の深さを計算することができる。具体的には、基板表面からの深さとその深さ位置における転位長さとを比較し、転位長さLd>深さDの関係を満たす場合、その関係を満たす深さにBMDが存在すれば、BMDから発生した転位が基板表面にまで伸びてスリップを発生させる可能性を示している。例えば、図19に示すエネルギー条件2における転位長さ曲線では、基板表面からの深さが5μmのとき、転位長さは約6.5μmであり、上記転位長さLd>深さDの関係を満たすので、基板表面からの深さ5μmの位置にBMDが存在すると、発生した転位が基板表面にまで到達する可能性が非常に高いことが判る。このように構成された本発明の転位予測方法では、上記各ステップを経ることにより、LSA工程でのシリコン基板に含まれるBMDから発生する転位の長さを正確に予測することができる。これによりBMDから発生する転位を最小限にすることが可能なLSA処理条件を探し出すことができる。また、BMDから発生する転位が基板表面に達しない範囲や基板表面に達する密度が小さい範囲のBMD核密度、酸素濃度を探し出すことができる。従って、このシミュレーションを用いて得られた予測結果に基づき、BMDから発生する転位を最小限にすることが可能な条件をLSA処理条件として選定することで、LSAによる転位の発生を抑制することができる。このシミュレーションを用いて得られた予測結果に基づき、BMDから発生する転位が基板表面に達しない範囲のBMD核密度、並びに酸素濃度を有するシリコン基板をLSA処理対象基板として選定することで、LSAによる転位の発生を抑制することができる。このシミュレーションを用いて得られた予測結果に基づき、BMDから発生する転位の基板表面に達する密度が小さい範囲の、BMD核密度、並びに酸素濃度を有するシリコン基板をLSA処理対象基板として選定することで、LSAによる転位の発生を抑制することができる。次に本発明の実施例を詳しく説明する。<実施例1>第1ステップ先ず、図4に示すようにシリコン基板の断面の一部分、幅X方向は4.5mm、高さY方向は775μmの構造を持つ計算範囲を複数のメッシュ要素に細分化してメッシュ化した。なお、高さY方向の計算領域は300mmシリコン基板の一般的な厚みに相当する。このときに、X軸の方向は[110]方向と、Y軸の方向は[100]方向とそれぞれ一致する面を仮定した。また、各メッシュ要素の幅は一定とし、また、表面部をより詳細に計算するために、表面に近づくほどメッシュ要素の高さを小さくした。表面部でのメッシュ要素最小サイズを高さ約1.2μm×横幅15μmとした。このサイズでは120μm幅のレーザービームが横幅15μmのメッシュ要素8個分と一致するようになる。メッシュ化されたシリコン基板の計算範囲における横方向のメッシュ要素数は300個、深さ方向のメッシュ要素数は100個、総メッシュ要素数は30000個である。第2ステップ次いで、レーザービームのスキャン速度、レーザービームのサイズ、レーザービームのエネルギー及びスキャン経路から、レーザービームが照射されるメッシュ要素のエネルギーの時間変化を計算した。レーザービーム条件として、レーザービームのサイズを幅120μm、スキャン速度を200mm/秒に設定した。なおこのサイズ条件から、上記第1ステップのシミュレーションモデル構造で説明した要素サイズ幅を決定している。また、レーザービームパワーは、表面熱流速のエネルギー条件12.4×104J/cm2s、エネルギー条件22.6×104J/cm2s及びエネルギー条件32.8×104J/cm2sの3種類を選定した。このエネルギーレベルは様々なテストから、LSA工程で実際に行われているアニール温度である1100~1300°Cになるような水準を選定している。また、実際のLSA工程時間に該当するスキャン時間は、一定な温度に大体に推移する時間を考慮して11ミリ秒とした。第3ステップ次に、各メッシュ要素の座標データ、物性値、基板の初期温度及び上記第2ステップで計算されたエネルギーの時間変化をそれぞれコンピュータに入力した。物性値は、シリコンの密度2.3g/cm3、ヤング率1.7GPa、ポアソン比0.26とし、熱伝導度、比熱、熱膨張係数は、図8~図10に示す温度依存性データを入力データに使用した。また、基板の初期温度は400°Cで、伝熱計算の初期条件として設定された。第4ステップ次に、基板の初期温度、熱伝導率に基づいたレーザービームの入射エネルギーを負荷として、基板の温度分布の時間変化を求めた。図11に、レーザーエネルギーとして与えたエネルギー条件12.4×104J/cm2s、エネルギー条件22.6×104J/cm2s、エネルギー条件32.8×104J/cm2sにおける基板の温度の時間変化を示す。図11から明らかなように、エネルギー条件12.4×104J/cm2sでは最終温度が約1140°Cまで、エネルギー条件22.6×104J/cm2sでは約1230°Cまで、エネルギー条件32.8×104J/cm2sでは約1300°Cまで達していることが判る。図12に、エネルギー条件22.6×104J/cm2s、スキャン時間11ミリ秒後におけるメッシュ化した基板の温度分布を示す。図12に示す温度分布では、基板表面のレーザービームが照射される部分が約1230°Cにまで上昇し、その上昇した高温域が深さ方向に拡散する模様を示している。第5ステップ次に、上記第4ステップで求められた温度分布を利用して応力計算を実施し、基板に生じる熱応力分布の時間変化を求めた。図13に、エネルギー条件22.6×104J/cm2s、スキャン時間11ミリ秒後におけるメッシュ化した基板の応力分布を示す。図13aに示すX方向の引張応力Tensilestress成分σxxでは、基板表面側では最大約241.7MPaの圧縮Compressive応力がかかっている。一方、基板裏面側には、圧縮応力に比べて小さい、最大約37.06MPaの引張Tensile応力がかかっている。図13bに示すX方向又はY方向せん断応力Shearstress成分τxyは、上記σxx成分の応力に比べて約0.1倍の大きさであることと、ビームの前方では圧縮応力が、ビームの後方では引張応力が掛かり、ビームが通過することによって反転することが判る。なお、ここではX、Y方向を等方性に仮定したのでX方向、Y方向は同じ結果となる。図13cに示すY方向の引張応力成分σyyでは、最大値が約21.5MPaのY方向引張り成分が断面方向に分布していることが判る。これらの結果から、100MPaを越えるような大きな応力の発生は、表面部に限定して存在していて、この発生した応力が深さ分布によって存在するBMDに到達すると転位に繋がる可能性を示している。第6ステップ次に、メッシュ化された基板の任意の位置Aを設定し、上記第4ステップで求められた温度分布及び上記第5ステップの基板に生じた熱応力分布から、設定した任意の位置Aにおける深さ方向の分布を各時間で選び出した。図14に、メッシュ化された基板に設定したある任意の位置Aを示す。また、図15に、エネルギー条件22.6×104J/cm2sにおける、任意位置Aの各時間ごとに選び出した深さ方向の温度分布を示す。第7ステップ次に、上記第6ステップで選び出した温度分布と基板中に含まれるBMDサイズ及び残存酸素濃度から、BMDから発生する転位の臨界応力を求めた。フィッティング式は図16に示す実験結果を基に導き出された上記式4に示す板状析出物に対する関連フィッティング式を用いた。また、残存酸素の濃度は1.1×1018/cm3、析出物サイズは100nmと仮定した。図17に、エネルギー条件22.6×104J/cm2sにおける、BMDから発生する転位の臨界応力CriticalStressτBMDを示す。図17から明らかなように、任意位置Aにおいて選び出した計算の領域は約6ミリ秒でレーザービームが通過するので、その時に表面側の臨界応力が低くなっていて深さ方向に徐々に増加するのが判る。臨界応力の時間変化からは、ビームが通過する時間には温度が上昇し、それに従って臨界応力も減少するが、しかし、ビームが通過する前、又は通過した後は温度も低下し、臨界応力が急激に増加するのが判る。これで、ビーム通過時に温度上昇で臨界応力が低くなるので、BMDに臨界応力以上の大きさの応力が掛かれば転位が発生する可能性を示している。第8ステップ次に、上記第6ステップで選び出した熱応力分布の応力成分σxx、σyy、τxyから、シリコン単結晶の主スリップ面である111面の[110]方向のせん断応力成分τxyを計算した。第9ステップ次に、上記第8ステップのせん断応力τxyと上記第7ステップの臨界応力τBMDとを各深さ位置ごとに比較し、比較結果から各深さ位置の転位長さを求めた。図18に、各深さ位置におけるせん断応力τxyと臨界応力τBMDとを重ね合わせた図を示す。第10ステップ次に、上記第9ステップで求めた転位長さから、スリップ発生の深さを計算した。図19に、本発明のシミュレーション方法により計算された、エネルギー条件1~エネルギー条件3における転位の長さを示す。図19から明らかなように、温度が高くなる条件であるほど、転位の長さが長く伸びることが判る。具体的には、エネルギー条件12.4×104J/cm2sでは基板表面から約3μmの深さにBMDが存在すれば、基板表面にまで転位が伸びてスリップが起きる可能性を示している。また、エネルギー条件22.6×104J/cm2sでは基板表面から約6μmの深さに、エネルギー条件32.8×104J/cm2sでは基板表面から約9μmの深さに、それぞれBMDが存在すれば、スリップに繋がり、オーバーレイ位置合わせ誤差の問題を起こす可能性があることを示している。"