(57)【要約】【課題】精度が高く、計算量が少なく、かつ確実に、結晶基板中の目的物質の拡散を予測することができる目的物質の拡散予測方法を提供する。【解決手段】本発明の目的物質の拡散予測方法は、アニールを行ったときの結晶基板中の目的物質の拡散を計算によって予測する方法であって、複数の拡散種を想定し、かつ、アニール時間を複数の時間領域に分け、前記各時間領域毎に、それぞれ、前記各拡散種を第１要素と第２要素とに振り分けるとともに、計算における時間刻み間隔Δｔを設定し、アニール温度Ｔに応じて、目的物質の拡散を予測するための時間に関するパラメーターを設定し（ステップＳ１０１〜Ｓ１０８）、これらの設定値に基づいて計算することを特徴とする。【選択図】図１１
2005057143.tif 000002
【特許請求の範囲】
【請求項１】
アニールを行ったときの結晶基板中の目的物質の拡散を計算によって予測する方法であって、
複数の拡散種を想定し、かつ、アニール時間を複数の時間領域に分け、前記各時間領域毎に、それぞれ、前記各拡散種を第１要素と第２要素とに振り分けるとともに、計算における時間刻み間隔Δｔを設定し、
アニール温度Ｔに応じて、目的物質の拡散を予測するための時間に関するパラメーターを設定し、これらの設定値に基づいて計算することを特徴とする目的物質の拡散予測方法。
【請求項２】
アニールを行ったときの結晶基板中の目的物質の拡散を計算によって予測する方法であって、
複数の拡散種を想定し、かつ、アニール時間を複数の時間領域に分け、前記各時間領域毎に、それぞれ、前記各拡散種を第１要素と第２要素とに振り分けるとともに、計算における時間刻み間隔Δｔを設定し、
アニール温度Ｔに応じて、前記各時間領域の時間帯を設定し、これらの設定値に基づいて計算することを特徴とする目的物質の拡散予測方法。
【請求項３】
アニールを行ったときの結晶基板中の目的物質の拡散を計算によって予測する方法であって、
複数の拡散種を想定し、かつ、アニール時間を複数の時間領域に分け、前記各時間領域毎に、それぞれ、前記各拡散種を第１要素と第２要素とに振り分けるとともに、前記各時間領域毎に、それぞれ、アニール温度Ｔに応じて、計算における時間刻み間隔Δｔを設定し、これらの設定値に基づいて計算することを特徴とする目的物質の拡散予測方法。
【請求項４】
アニール温度Ｔに応じて、前記各時間領域の時間帯を設定する請求項３に記載の目的物質の拡散予測方法。
【請求項５】
前記複数の時間領域には、前記時間刻み間隔Δｔが一定に設定される時間領域と、前記時間刻み間隔Δｔがアニール経過時間ｔに応じて増大するように設定される時間領域とが含まれる請求項１ないし４のいずれかに記載の目的物質の拡散予測方法。
【請求項６】
前記第１要素に振り分けられる拡散種は、平衡状態に達していないものとして扱われる拡散種であり、前記第２要素に振り分けられる拡散種は、平衡状態に達したものとして扱われる拡散種である請求項１ないし５のいずれかに記載の目的物質の拡散予測方法。
【請求項７】
前記第１要素に振り分けられた拡散種の濃度は、第１の関係式を用いて算出し、前記第２要素に振り分けられた拡散種の濃度は、第２の関係式と、前記算出した前記第１要素に振り分けられた拡散種の濃度とを用いて算出する請求項１ないし６のいずれかに記載の目的物質の拡散予測方法。
【請求項８】
前記第１の関係式は、下記式１である請求項７に記載の目的物質の拡散予測方法。
【数１】
2005057143.tif 000009
【請求項９】
前記第２の関係式は、下記式２である請求項７または８に記載の目的物質の拡散予測方法。
【数２】
2005057143.tif 000010
【請求項１０】
前記目的物質は、ホウ素である請求項１ないし９のいずれかに記載の目的物質の拡散予測方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、目的物質の拡散予測方法に関するものである。
【背景技術】
【０００２】
従来、ホウ素等の目的物質の結晶基板中への拡散を予測するシミュレータが知られている（例えば、特許文献１参照）。
半導体デバイスは、小型化によって高性能化が達成されており、今後ますます微細化されると予測できる。微細プロセス設計のための物性予測の精度向上が求められる一方で、０．１μｍ設計ルール以降は物性が急変し、予測が困難になる。結晶基板中の目的物質の拡散、特にホウ素拡散に関しては、より高濃度に、極浅の領域にホウ素を留まらせることが要求される。
しかしながら、従来のシミュレータでは、濃度予測精度が悪い。また、精度向上のため数多くの実験を行う必要が生じ、例えば、トランジスタ、メモリー、ＬＳＩ、ＣＰＵ等の半導体デバイスの設計開発には、長期間、大規模の投資が必要になる。
【０００３】
【特許文献１】特開平５−１１４５６９号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
本発明の目的は、精度が高く、計算量が少なく、かつ確実に、結晶基板中の目的物質の拡散を予測することができる目的物質の拡散予測方法を提供することにある。
【課題を解決するための手段】
【０００５】
このような目的は、下記の本発明により達成される。
本発明の目的物質の拡散予測方法は、アニールを行ったときの結晶基板中の目的物質の拡散を計算によって予測する方法であって、
複数の拡散種を想定し、かつ、アニール時間を複数の時間領域に分け、前記各時間領域毎に、それぞれ、前記各拡散種を第１要素と第２要素とに振り分けるとともに、計算における時間刻み間隔Δｔを設定し、
アニール温度Ｔに応じて、目的物質の拡散を予測するための時間に関するパラメーターを設定し、これらの設定値に基づいて計算することを特徴とする。
これにより、任意のアニール温度Ｔ、特に、半導体のプロセスに使用される温度領域（例えば、０〜１２００℃程度）のすべてにおいて、計算が破綻することなく、精度良く、迅速かつ確実に、目的物質の濃度分布を予測することができる。
【０００６】
本発明の目的物質の拡散予測方法は、アニールを行ったときの結晶基板中の目的物質の拡散を計算によって予測する方法であって、
複数の拡散種を想定し、かつ、アニール時間を複数の時間領域に分け、前記各時間領域毎に、それぞれ、前記各拡散種を第１要素と第２要素とに振り分けるとともに、計算における時間刻み間隔Δｔを設定し、
アニール温度Ｔに応じて、前記各時間領域の時間帯を設定し、これらの設定値に基づいて計算することを特徴とする。
これにより、任意のアニール温度Ｔ、特に、半導体のプロセスに使用される温度領域（例えば、０〜１２００℃程度）のすべてにおいて、計算が破綻することなく、精度良く、迅速かつ確実に、目的物質の濃度分布を予測することができる。
【０００７】
本発明の目的物質の拡散予測方法は、アニールを行ったときの結晶基板中の目的物質の拡散を計算によって予測する方法であって、
複数の拡散種を想定し、かつ、アニール時間を複数の時間領域に分け、前記各時間領域毎に、それぞれ、前記各拡散種を第１要素と第２要素とに振り分けるとともに、前記各時間領域毎に、それぞれ、アニール温度Ｔに応じて、計算における時間刻み間隔Δｔを設定し、これらの設定値に基づいて計算することを特徴とする。
これにより、任意のアニール温度Ｔ、特に、半導体のプロセスに使用される温度領域（例えば、０〜１２００℃程度）のすべてにおいて、計算が破綻することなく、精度良く、迅速かつ確実に、目的物質の濃度分布を予測することができる。
【０００８】
本発明の目的物質の拡散予測方法では、アニール温度Ｔに応じて、前記各時間領域の時間帯を設定するのが好ましい。
本発明の目的物質の拡散予測方法では、前記複数の時間領域には、前記時間刻み間隔Δｔが一定に設定される時間領域と、前記時間刻み間隔Δｔがアニール経過時間ｔに応じて増大するように設定される時間領域とが含まれるのが好ましい。
【０００９】
本発明の目的物質の拡散予測方法では、前記第１要素に振り分けられる拡散種は、平衡状態に達していないものとして扱われる拡散種であり、前記第２要素に振り分けられる拡散種は、平衡状態に達したものとして扱われる拡散種であるのが好ましい。
本発明の目的物質の拡散予測方法では、前記第１要素に振り分けられた拡散種の濃度は、第１の関係式を用いて算出し、前記第２要素に振り分けられた拡散種の濃度は、第２の関係式と、前記算出した前記第１要素に振り分けられた拡散種の濃度とを用いて算出するのが好ましい。
【００１０】
本発明の目的物質の拡散予測方法では、前記第１の関係式は、下記式１であるのが好ましい。
【００１１】
【数１】
2005057143.tif 000003
【００１２】
本発明の目的物質の拡散予測方法では、前記第２の関係式は、下記式２であるのが好ましい。
【００１３】
【数２】
2005057143.tif 000004
本発明の目的物質の拡散予測方法では、前記目的物質は、ホウ素であるのが好ましい。
【発明を実施するための最良の形態】
【００１４】
以下、本発明の目的物質の拡散予測方法（拡散シミュレーション方法）を添付図面に示す好適実施形態に基づいて詳細に説明する。
本発明の目的物質の拡散予測方法は、アニールを行ったときの結晶基板中の目的物質の拡散を計算によって予測する方法、すなわち、結晶基板中の目的物質の濃度分布を計算によって予測する方法である。この本発明の理解を容易にするため、まずは、本発明の前提や基礎知識について説明し、その後、本発明の要部（特徴）について説明する。なお、以下の説明では、前記本発明の前提や基礎知識についての説明においても、それを、本発明（本発明の一部）として説明する。
【００１５】
図１は、計算機シミュレーションを用いたプロセス設計の一般的アプローチを示す図である。量子力学計算を使い、第１原理から、所望の各種プロセス条件について構成される可能性のある材料の電子的状態と構造を予測する。
電子機器プロセスのマルチスケールのモデル化方法２００を図２に示す。第１のステップにおいて、プロセス中に起こるかもしれない一つまたは複数の微細反応過程が特定される（ステップ２１０）。反応過程には、シリコン（Ｓｉ）格子へのホウ素（Ｂ）注入等、電子材料のプロセスにおいて使用される従来の技術を含めることができる。微細反応過程には、クラスターや複合体を形成するドーパント、基板原子、不純物の会合、これらの複合体、ドーパント要素の拡散、複合体の解離等を含め、作業中に発生する、または発生するかもしれない反応過程が含まれる。次に、所望の微細反応過程に関する基本データを入手する（ステップ２２０）。このステップでは、量子力学（ＱＭ）を使って、各種状態のエネルギー過程と構造、原子と欠陥が部位間を遷移する各種反応過程に関する障壁、拡散、会合、解離速度の予測に必要なその他のデータを予測することができる。これは、各種の定常構造の相対的エネルギーと各部位を接続する力学的、動力学的反応過程についての活性化エネルギーとエントロピーを計算することによって行われる。たとえば、活性化エントロピーは、結合のための状態および遷移状態の振動密度から得ることができる。さらに、バンド状態と電気的動作のＱＭ計算により、ドーパント、欠陥、不純物のさまざまな置換、格子間介入、クラスターに関する電子的特性（たとえば、導電性、ホール効果、磁気抵抗、熱電能）の予測に関連するパラメータが得られる。量子力学計算は、たとえば、既知のＣＰＭＤ  Ｖ３．３パッケージ(UJ. Hutter, A. Alavi, T. Deutch, M. Bernasconi, S. Goedecker, D. Marx, M. Tuckerman, M. Parrinello, MPI fuer FestkoerperforschungおよびIBM Zurich Research Laboratory 1995-1999)で実現されている密度汎関数理論、あるいはその他周知の技術（たとえば、Accelrys Inc. (サンディエゴ)が販売するCastepおよびＤＭＯＬ３、オレゴン州ポートランドのSchrodinger Inc.が販売するJaguar等のプログラム）を使用して実行できる。計算の初期パラメータの仕様は、モデル化する系によって異なる。ある場合には、既知の結晶構造からスタートし、統計的手法（たとえば、モンテカルロサンプリング）を使って欠陥と不純物の分布を代入することによって生成できる。別の場合には、加熱冷却プロセスに力場を用いた分子動力学（ＭＤ）を使って、系のランダム化とアニールを行うことができる。これらの手法は、構造の組み合わせを生成するのに使用でき、この組み合わせはＱＭとＭＤを使って計算し、その後の分析のための重要な構造を決定することができる。
【００１６】
ＱＭ計算またはＭＤ計算あるいはステップ２２０の別の場所から得られた基本反応レートパラメータデータは、分子動力学を使って、数種の、または多種の欠陥と不純物分布を取り入れた集合系の特性をモデル化するのに使用できる（ステップ２３０）。このステップにおいて、実用的に、または実験的に興味がある系を模倣するための、系のアンサンブルを構成する目的で、欠陥や不純物を適当な方法で（例えば、本来は実験的に設定される条件を模倣する方法）空間的に分布させたたくさんの系を考えることができる。これらの系は、前述のように構成できる。次に、これらの計算はＱＭその他からの基本反応レートパラメータデータを使い、上記の欠陥と不純物の拡散を含むミクロ、メソ、マクロ系の反応レートを予測し、欠陥とドーパントの力学を判断する。
【００１７】
ある実施形態において、これらの計算は、キネティックモンテカルロ（ｋＭＣ）動力学を使って実行される。この計算方法では、複数の時間スケールを許して独立した進度で系を時間発展させることで、多様なプロセス条件を記述できる実用的な手段を提供している。マルチスケールｋＭＣシミュレーションは、欠陥／ドーパント反応の長時間経過後の結果に関する情報を提供し、どのクラスターが主要な要素を構成するかを示し、重要な反応経路を特定する。
【００１８】
基本データ（たとえば、ステップ２２０のＱＭ計算から得る）と反応情報（たとえば、ステップ２３０のｋＭＣ力学計算から得る）を使い、ドーパントと欠陥要素の分布とクラスターを予測する（ステップ２４０）。これらの要素の平衡および非平衡分布は、初期濃度、温度およびその他外的条件（たとえば、応力、電圧、磁場、電磁界、温度変化、外的放射）に関して予測される。予測された分布は、時間とプロセス条件に関して、上記の原子的アプローチ（ＱＭ、ＭＤ）のほかに、原子の状況に基づき連続体近似をした偏微分方程式を有限要素法（ＦＥＭ）がかかわる階層的な方法で解かれる。
【００１９】
プロセスには、上述のようなドーパント、欠陥、クラスターまたは複合体を含む、多くの要素がかかわる可能性があり、その濃度は、所与の時間において相互依存することがある。したがって、たとえば、シリコン格子にホウ素のイオン注入を行うことで、置換または格子間ホウ素原子、格子間シリコン原子、空孔、さまざまな数のホウ素およびシリコン原子のクラスターを含む多くの種が作られる。これらの種の一部は、ここでは「高速」または「超高速」種または要素と呼ぶものとし、その拡散係数または解離速度（それぞれ拡散エネルギー障壁と結合エネルギーによって予測可能）が非常に高く、これらの種の局所濃度は、これらが他の要素から遠くに拡散する、または他の要素に解離するにつれ、時間とともに大きく異なる。
【００２０】
系が時間発展すると、一部の高速プロセスによって、相対的濃度がほぼ一定で、全体的な母集団が時間とともにゆっくりと変化する定常状態が実現する。このような状況では、原子の動力学は、原子的な方程式から判断されるパラメータの構成式に置き換えることができる。これによって、より長時間の現象をメソスケールまたは連続体方程式で表すことができ、実践的な時間と長さスケールを大きく伸ばすことができる。
【００２１】
このように、高速種の濃度は次のように表すことができる。
∂Ｃｋ／∂ｔ＝∇（Ｄｋ∇Ｃｋ）＋ＫＲ，ｉｊＣｉＣｊ＋ＫＤ，ｍＣｍ−ＫＲ，ｋｌＣｋＣｌ−ＫＤ，ｋＣｋ
ただし、ＤｋとＣｋはそれぞれ、要素ｋの拡散係数と濃度であり、反応および分解係数は次式で表される。
【００２２】
ＫＲ，ｉｊ＝４πＲＰ（Ｄｉ＋Ｄｊ）
ＫＤ，ｋ＝ｎＨ／４ＫＲ，ｉｊｅｘｐ（−Ｅｂ，ｋ／ｋＢＴ）
ただし、Ｒｐは捕獲半径、ｎＨは格子位置の濃度、Ｅｂはｋからｉ，ｊへの解離エネルギーである。上式において、右辺第１項は拡散、第２項、第３項はそれぞれ、ｉ＋ｊクラスター化（ｉ＋ｊ→ｋ）とｍ分解（ｍ→ｋ＋ｌ）によって生成される。ｋ要素は、第４項（ｋ→ｉ＋ｊ）と第５項（ｋ＋ｌ→ｍ）で示されるような、逆プロセスによって取り除かれる。
【００２３】
アニーリング（アニール）等のプロセス開始時に、高速要素は、他の低速で変化する種（「低速」要素。すなわち結合エネルギーが高く、拡散係数が小さい種）との熱平衡から離れている。これらの高速要素の動力学を説明するには、上記の偏微分方程式を非常に小さな時間刻みで解かなければならない（最も短いシミュレーション時間スケールより小さいか少なくとも同程度で、最も短いシミュレーション時間スケールの近似値はΔｘ２／Ｄまたは１／ｋｄで得られる。ただし、Δｘはシミュレーションの空間離散化サイズとする）。したがって、たとえば、１０００℃の時、シリコン格子中の格子間単一Ｓｉ、Ｂｅ−Ｓｉｉ対または空孔の動力学を説明するのに、１０−１２秒のオーダーでの時間刻みが必要となる。上記のように小さな時間刻みを使って要素の濃度が減少することを説明することは、１秒より長い代表的なシミュレーション時間スケールにとって計算集約的で、多くの高速要素がかかわる一般的な代表的シミュレーションに、何時間、あるいはそれ以上の時間がかかる。
【００２４】
しかしながら、高速要素は通常、低速要素との熱平衡、それらの濃度が基本的に低速要素のそれによって決定される、いわゆる擬似的定常状態（平衡状態）に達すると認識することによって、計算時間は大幅に短縮できる。この段階で、高速要素の濃度は、以下の偏微分方程式を解くことによって決定できる。
∇（Ｄｋ∇Ｃｋ）＋ＫＲ，ｉｊＣｉＣｊ＋ＫＤ，ｍＣｍ−ＫＲ，ｋｌＣｋＣｌ−ＫＤ，ｋＣｋ＝０
つまり、所与の高速要素について擬似的定常状態に達する前の時点で、その高速要素の濃度変化は、前述のような小さな時間刻みで偏微分方程式を解くことによって示すことができ、高速要素が擬似的定常状態に達した後は、その濃度は、低速で変化する要素がかかわる偏微分方程式を解くことによって決定できる。
【００２５】
この段階で、時間刻みは次に最速の要素に合うよう（たとえば、次に最速の要素の拡散係数と解離速度による）、調整することができる。擬似的定常状態に達する高速要素を分離することで、時間刻みを連続的に増加できる。これに対し、低速要素の濃度は、適正な時間刻みで関連する偏微分方程式を単純に解くことによって決定され、これは低速要素の大きな結合エネルギーと小さい拡散係数のためにずっと大きく、したがって高速要素の場合ほど計算集約的とならない。
【００２６】
このように、計算時間は大幅に短縮することができ、その結果、ピコ秒から分までの時間スケールがかかわる非平衡プロセスを計算効率のよい（費用対効果の高い）方法で判断、予測することができる。さらに、これらの計算効率により、実際の装置のプロセス条件を最適化する上で不可欠な、大規模な二次元および三次元シミュレーションが可能となる。
【００２７】
[例−超極浅接合プロセス]
上述のようなマルチスケールシミュレーション方法を超極浅接合プロセス、特にシリコンへのホウ素の超極浅蒸着への応用について、図３〜図９を参照しながら説明する。この例では、ある一つの部位における複数の原子または欠陥のクラスターを含め、各種の置換または格子間部位におけるホウ素、シリコン、空孔の拡散を検討する。超極浅接合プロセスのためのマルチスケールモデル化の適用を図３に示す。
【００２８】
[１．微細反応過程：量子力学]
まず、量子力学計算３１０を使い、関連する微細反応過程をモデル化した。すべての原子構造は、ＣＰＭＤ  Ｖ３．３パッケージにおいて用いられているように、局所密度近似（ＬＤＡ）（D.M. CeperleyとB.J. Alder, Phys. Rev. Lett. 45, 566 (1980), J. PerdewとA. Zunger, Phys. Rev. B 23, 5048 (1981)）を使い、密度汎関数理論（ＤＦＴ）により最適化された。一般化勾配近似（ＧＧＡ）（J.P. Perdew et al., Phys. Rev. Lett. 77, 3865 (1996)）を使ったＬＤＡ計算で得られた構造のエネルギー過程の計算でも同じ結果が得られた。なお、非局所的、ノルム保存型擬ポテンシャル（S. Goedecker et al., Phys. Rev. B 54, 1703 (1996)）と２０Ｒｙの平面波カットオフエネルギーを使った。検討対象とした欠陥系は、純粋シリコン中に２．３５ÅのＳｉ−Ｓｉ結合距離を実現する一定体積での１９２原子を含む超格子を使ってモデル化された。すべての原子位置は、すべての残留力が５×１０−４Hartree/Bohrより小さくなるまで、十分に緩和された。大きなスーパーセルであることから、ブリルアンゾーン（ＢＺ）積分にガンマ点（ｋ＝０の点）のみを使う。もちろん、別の系では詳細（たとえばｋ点の数）を変えることもできる。一般に、プロトタイプ系についてのこれらのパラメータの範囲は、容認できるコストで容認できる精度を提供できるパラメータとされる。前述のとおり、多くの計算プログラムを使って、詳細なＱＭ計算を行うことができる。
【００２９】
まず、これらの方法（ＱＭ、ＭＤ等）を使い、接合プロセスの総合的なモデルを作る上で重要な、欠陥と欠陥、欠陥とドーパントの相互作用をより十分に理解した。特に、ここで検討している例については、シリコン内の超極浅ｐ型接合のためのホウ素ドーピングによって、ホウ素ホウ素およびホウ素シリコン格子間クラスター化と関連する多くの複雑な拡散現象がどのようにもたらされるかを理解する必要があった。超極浅接合製造には非常に高いドーピング濃度が必要となるため（たとえば、１０２０ｃｍ−３以上）、このようなクラスター化は不可避であろう。シリコンクラスターはホウ素の過渡的増速拡散（ＴＥＤ）の主原因となる単一シリコンの格子間原子の発生源として知られているため、イオン注入中およびアニーリング開始時のシリコン格子間クラスターの形成に関する具体的モデルを作ることが重要である。先に、量子力学計算（密度汎関数理論）で主要基本プロセスを特定できると説明した。このような研究は、本願の手法の最重要部分である。特に、本発明による系の例においては、次の（１）〜（３）の基本プロセスを決定した。
（１）ホウ素ホウ素対（置換−格子間、Ｂｓ−Ｂｉ）の拡散
（２）Ｂｓ−Ｂｓ−Ｂｉ複合体の分解
（３）シリコン自己格子間クラスター化におけるホウ素の触媒としての役割
【００３０】
[Ｂｓ−Ｂｉ拡散]
ＱＭを使い、ホウ素対の拡散が、超極浅接合プロセス中の拡散プロファイルを確定する上で重要な役割を果たすことがわかった（高濃度のホウ素と高温のアニーリングが必要）。図４は、ホウ素対が、［Ｂ−Ｂ］ｓ−＜１１１＞，Ｂｓ−Ｂｉ，Ｂｓ−Ｂｓ−Ｓｉｉ（Ｓｉｉは格子間シリコン）と表される３つの局所における最低状態を通じて、［Ｂ−Ｂ］ｓ−＜００１＞の一つの最低エネルギー構造から隣接する等価部位での等価構造へと拡散する反応経路を示す。このように、切れたり、歪んだりしている結合の数をＢ−Ｂ拡散中、最低に保つことができる。拡散の活性化エネルギーは、局所密度近似で１．８７ｅＶと予測される（一般化勾配近似では１．８１ｅＶ）。
【００３１】
[Ｂｓ−Ｂｓ−Ｂｉ分解]
Ｂｓ−Ｂｓ−Ｂｉ複合体は、シリコン内にクラスターを含むホウ素の中で最も安定していることが知られている。Ｂｓ−Ｂｓ−Ｂｉ複合体は、Ｂｓ−ＢｓとＢｉに分解できる。つまり、可動ホウ素原子は逃げ、Ｂｓ−Ｂｓが背後に残り、Ｂｓ−Ｂｓ−ＢｉをＢｓ−Ｂｓ＋Ｂｉに解離するエネルギーは約４．０ｅＶとなる。また、Ｂｓが最も近い部位に拡散した結果として、Ｂｓ−Ｂｓ−ＢｉはＢｓ−ＢｉとＢｓに解離できる。図５は、Ｂｓ−Ｂｓ−Ｂｉ拡散のメカニズムとして考えられるものを示す（Ｂｉ拡散は（ａ）、Ｂｓ拡散は（ｂ））。さらに、計算は、格子間Ｓｉを捕獲してＢｓ−Ｂｓ−ＢｉからＢｓ−Ｂｉ−Ｂｉに変化することは、３つのホウ素からなるクラスターを分解するよりずっと低い熱消費となり、Ｂｓ−Ｂｉ−ＢｉをＢｓ−Ｂｉ＋Ｂｉに解離する１．０ｅＶのエネルギーは、Ｂｓ−Ｂｓ−Ｂｉに関する４．０ｅＶの場合よりずっと低いことを示唆している。これは、Ｓｉ格子間原子の数の増加により、ホウ素クラスターの分解が促進された実験結果と一致している。
【００３２】
[自己格子間クラスター化におけるホウ素の触媒としての役割]
一つの格子間Ｓｉ原子は、格子置換ホウ素により捕獲されてＢｓ−Ｓｉｉ対を形成する。Ｂｓ−Ｓｉｉの結合エネルギーを約１．０ｅＶと算出する。Ｂｓ−Ｓｉｉ対はさらに、別の単一Ｓｉ格子間原子を捕獲することによって安定化される。Ｂｓ−ＳｉｉとのＳｉｉの結合エネルギーを約１．６ｅＶと算出する。このようなＢ−Ｓｉｉ，ｘ複合体は、可動Ｓｉ格子間原子を捕獲、安定化する。その一方で、Ｂの結合エネルギーは、クラスター内部の格子間Ｓｉ原子の数が増すにつれて弱まる。これは、ホウ素原子が事実上、低い結合エネルギー障壁を乗り越えてクラスターから逃れることができることを示している。全体として、ホウ素原子はＳｉ格子間クラスター化の触媒としての役割を果たす。図６は、電気的中性のＢ−Ｓｉｉ，ｘ複合体（ｘ＝１−４、（ａ）Ｂｓ−Ｓｉｉ、（ｂ）Ｂｓ−Ｓｉｉ，２、（ｃ）Ｂｓ−Ｓｉｉ，３、（ｄ）Ｂｓ−Ｓｉｉ，４、（ｅ）Ｓｉｉ，４、（ｆ）Ｓｉｉ，３）と、単一Ｓｉ格子間原子とホウ素の結合エネルギーを示す図である。この所見から、アニーリングの初期段階で小さなシリコンクラスターが形成される様子を知ることができる。
【００３３】
[２．欠陥−ドーパントの力学：キネティックモンテカルロ]
再び図３を参照すると、キネティックモンテカルロ分子動力学シミュレーション３２０を使い、ドーパントの拡散と欠陥−ドーパント相互作用が説明されている。これらの計算において、欠陥、ドーパント、クラスターの位置は、１０−１２秒から１０２秒までのマルチスケールシミュレーションの時間スケールの時間に関連して追跡することができ、１台のプロセッサについてのすべて包含したシミュレーションが数分で行われる。構成要素はすべて基本的な拡散係数と解離速度の基本的属性を有する質点として扱われる。主としてab initio ＱＭ計算に基づき、これらのｋＭＣシミュレーションに関するデータ（たとえば、可動種の温度依存拡散係数とクラスターの結合エネルギー）を得た。拡散係数の温度依存は次のように表すことができる。
【００３４】
Ｄｄｉｓｓ＝Ｄｏｅｘｐ（−Ｅｍ／ｋＢＴ）
ここで、Ｄｄｉｆｆは拡散係数、Ｄｏはプレファクタ、Ｅｍは拡散エネルギー障壁、ｋＢはボルツマン定数、Ｔは基板温度である。同様に、解離速度は次式で合わされる。
Ｄｄｉｓｓ＝Ｄｏｅｘｐ（−（Ｅｍ＋Ｅｂ）／ｋＢＴ）
ただし、Ｅｂは離れる可動種と残りの構成要素との結合エネルギーである。
【００３５】
図７は、ｋＭＣシミュレーションの略図である。このシミュレーションにおいて、クラスター化と解離を含む各種の欠陥−ドーパントの動力学が起こる。まず、各時間刻みで、シミュレーションボックスの中で起こりうるすべての事象を選択する（ステップ７１０）。次に、各粒子の状態を更新する（ステップ７２０）。これらのステップは、シミュレーション時間の合計が所定の量になるまで（ステップ７３０の「イエス」）繰り返す（ステップ７３０の「ノー」）。
【００３６】
図８は、ｋＭＣシミュレーションの結果を示す。１０１９ｃｍ−３の濃度でＢｓ−Ｓｉｉ対から始め、９００ＫでのｋＭＣシミュレーションを実行した。１μ秒後、Ｂｓ−Ｓｉｉの多くは解離し、Ｓｉ格子間原子がクラスターの形成を始めた。シリコンのクラスターは大きくなり、アニール時間でより安定した状態となった。最終的に小さなクラスターはすべて消失して大きなものが形成され、オストワルト成長が見られた。このｋＭＣシミュレーションは、高温処理における欠陥とドーパントの挙動の様子を示す。原子スケールのプロセスをこのように視覚化することは、欠陥−ドーパントの動力学を理解するのに役立つ。しかしながら、このようなｋＭＣアプローチには、１００から１０２１ｃｍ−３の範囲の種の濃度変化を考慮に入れながら、ドーピングプロファイルの時間発展を正確に説明する上で限界がある。マルチスケールシミュレーションにおけるこの問題を解決するために、以下に説明する連続体近似に基づくモデルも開発した。
【００３７】
[３．ドーピングプロファイルの時間発展：連続体力学]
図３に示すように、連続体力学３３０を使って、巨視的ドーピングプロファイルの時間発展をモデル化した。前述のように、一連の非平衡反応は次のように表すことができる。
∂Ｃｋ／∂ｔ＝∇（Ｄｋ∇Ｃｋ）＋ＫＲ，ｉｊＣｉＣｊ＋ＫＤ，ｍＣｍ−ＫＲ，ｋｌＣｋＣｌ−ＫＤ，ｋＣｋ
ただし、ＤｋとＣｋはそれぞれ、要素ｋの拡散係数と濃度であり、反応および分解係数は次のように得られる。
ＫＲ，ｉｊ＝４πＲＰ（Ｄｉ＋Ｄｊ）
ＫＤ，ｋ＝（ｎＨ／４）ＫＲ，ｉｊｅｘｐ（−Ｅｂ，ｋ／ｋＢＴ）
ただし、Ｒｐは捕獲半径、ｎＨはＳｉ格子部位の濃度（−５×１０２２ｃｍ−３）、Ｅｂはｋのｉとｊへの解離エネルギーである。
【００３８】
これらのシミュレーションにおいて、アニーリング開始時に、ホウ素はＢｓ−ＳｉまたはＢｓ−Ｓｉ２の形態で存在する。捕獲半径（Ｒｐ）は、２．５×１０−８ｃｍに設定される。自己格子間原子、Ｂ−Ｓｉｉ対とＢｓ−Ｂｉ対はそれぞれ、０．００５×ｅｘｐ（−０．４ｅＶ／ｋＢＴ）、０．００１×ｅｘｐ（−０．６８ｅＶ／ｋＢＴ）、１．５×ｅｘｐ（−２．４５ｅＶ／ｋＢＴ）ｃｍ２／秒で拡散できる。主要な欠陥−ドーパント複合体、Ｂｓ−Ｓｉｉ、Ｂｓ−Ｂｉ、Ｂｓ＝Ｂｓ−Ｂｉの拡散レートはそれぞれ、１０１３ｅｘｐ（−１．１５ｅＶ／ｋＢＴ）、３×１０１３ｅｘｐ（−２．１８ｅＶ／ｋＢＴ）、１０１７ｅｘｐ（−４．４６ｅＶ／ｋＢＴ）秒−１と予測される。その他のホウ素格子間シリコン複合体（シミュレーションで検討）にはＢｓ−Ｓｉｉ，２、Ｂｓ−Ｓｉｉ，３、Ｂｓ−Ｂｉ，２、Ｂｓ−Ｂｉ，３、Ｂｓ−Ｂｉ−Ｓｉｉ、Ｂｓ−Ｂｉ−Ｓｉｉ，２、Ｂｓ，２−Ｂｉ、Ｂｓ，２−Ｂｉ，２がある。
【００３９】
上述のように、超高速微細反応過程（拡散、クラスター化、解離）を考慮に入れて、マクロ的ドーピングプロファイルを説明するために使用されるアルゴリズムは、濃度が時間によって急速に変化する欠陥とドーパント要素に関する擬似的定常状態の概念を取り入れている。高温アニーリング開始時に、Ｓｉ自己格子間原子、空孔および、高い拡散係数または高い解離速度を有するＢｓ−Ｓｉｉ対をはじめとする多数の（平衡をはるかに超える）超高速要素がある。これらの要素の濃度が減少することを説明するには、非常に小さな時間刻みが必要となる。たとえば、１０００℃で、単一Ｓｉ格子間原子と空孔の動力学を説明するには、通常、１０−１２秒の時間刻みが必要である。また、Ｂｓ−Ｓｉｉ対の高い解離速度（〜１０９秒−１）により、Ｂｓ−Ｓｉｉ濃度変化を計算するのに、小さな時間刻み（１０−９秒よりはるかに短い）が必要となる。アニール温度が高くなると、時間刻みはさらに短縮しなければならない。しかし、１０−９から（最大）１０−６秒の間で、いわゆる高速要素は、その濃度が他の低速で変化する要素によって決定される擬似的定常状態に達する。つまり、擬似的定常状態に達する前は、高速要素の濃度変化は、小さな時間刻み（好ましくは、Δｘ２／Ｄｉおよび１／ｋｄ，ｉ、ただし、Δｘは空間離散化サイズ、Ｄｉとｋｄ，ｉは高速要素の拡散係数と解離速度）で偏微分方程式を解くことによって説明される。
【００４０】
例えば、Ｂｓ−Ｓｉについて、
ｄＣＢｓ−Ｓｉ／ｄｔ＝−ＤＢｓ−ｓｉｄ２ＣＢｓ−Ｓｉ／ｄｘ＋ｇＢｓ−ｓｉ−ｌＢｓ−ｓｉ
ただし、ｇＢｓ−ｓｉ＝ＤｓｉＲｐＣｓｉＣＢｓ＋ｋｄ，Ｂｓ−Ｓｉ２ＣＢｓ−Ｓｉ２＋ｋｄ，Ｂｓ−ＢｉＣＢｓ−Ｂｉ＋ｋｄ，Ｂｓ，２−Ｂｉ２ＣＢｓ，２−Ｂｉ，２
ｌＢｓ−ｓｉ＝ＤｓｉＲｐＣｓｉＣＢｓ−Ｓｉ＋ＤＢｓ−ＳｉＲｐＣＢｓ−Ｓｉ（ＣＢｓ＋ＣＢｓ−Ｂｉ＋ＣＢ−Ｂｉ２＋ＣＢｓｉ＋ＣＢｓ−Ｓｉｓ＋ＣＢｓ，２−Ｂ）＋ｋｄ，Ｂｓ−ＳｉＣＢｓ−Ｓｉ
ここで、Ｃｉ，Ｄｉ，ｋｄ，ｉは種ｉの濃度、拡散係数、解離速度である。しかしながら、高速要素が擬似的定常状態に到達すると、その濃度は低速で変化する要素にかかわる常微分方程式を解くことによって決定される。
【００４１】
例えば、Ｂｓ−Ｓｉについて、
０＝−ＤＢｓ−ｓｉｄ２ＣＢｓ−Ｓｉ／ｄｘ＋ｇＢｓ−ｓｉ−ｌＢｓ−ｓｉ
この段階で、時間刻みは、次に最速の要素に適合するよう（たとえば、次に最速の要素の拡散係数と解離速度による）調整することができる。擬似的定常状態に到達する高速要素を分離することにより、時間刻みは、継続的に増大できる。
【００４２】
次に、本発明の目的物質の拡散予測方法の要部（特徴）について説明する。なお、以下の説明では、結晶基板としてシリコン基板（単結晶シリコン基板）、目的物質としてホウ素を選択した場合、すなわち、シリコン基板中のホウ素の拡散の予測について、実施形態として説明する。
本発明の目的物質の拡散予測方法は、アニールを行ったときのシリコン基板（結晶基板）中のホウ素（目的物質）の拡散（濃度分布）を計算によって予測する方法である。
【００４３】
この場合、クラスター等の拡散種を複数想定する。本実施形態で、使用するクラスターは、全部で、BsI、BsI2、BsI3、BsB、BsB2、BsB3、BsBI、BsBI2、Bs2B、BsB2の１０種であり、これらのクラスターと、ホウ素とを拡散種として同時に扱う。
これにより、精度が高くなり、シリコン基板中のホウ素の拡散（濃度分布）を正確に予測する（求める）ことができる。
【００４４】
また、アニール時間を複数の時間領域に分ける。本実施形態では、アニール時間を４つの時間領域、すなわち、時間ｔ１、ｔ２、ｔ３の大小関係を、ｔ１＜ｔ２＜ｔ３としたとき、ｔ１以下（ｔ≦ｔ１）の第１の時間領域と、ｔ１を超え、ｔ２以下（ｔ１＜ｔ≦ｔ２）の第２の時間領域と、ｔ２を超え、ｔ３以下（ｔ２＜ｔ≦ｔ３）の第３の時間領域と、ｔ３を超える（ｔ＞ｔ３）第４の時間領域とに分ける。
【００４５】
また、各時間領域毎に、それぞれ、各拡散種を第１要素と第２要素とに振り分ける。本実施形態では、下記の通り振り分ける。
第１の時間領域（ｔ≦ｔ１）
第１要素：BsI、BsI2、BsI3、BsB、BsB2、BsB3、BsBI、BsBI2、Bs2B、BsB2
第１要素：なし
第２の時間領域（ｔ１＜ｔ≦ｔ２）
第１要素：BsI2、BsI3、BsB、BsB2、BsB3、BsBI、BsBI2、Bs2B、BsB2
第２要素：BsI
第３の時間領域（ｔ２＜ｔ≦ｔ３）
第１要素：BsB、BsB2、BsB3、BsBI、BsBI2、Bs2B、BsB2
第２要素：BsI、BsI2、BsI3
第４の時間領域（ｔ＞ｔ３）
第１要素：BsB、BsB2、BsB3、BsBI、BsBI2、Bs2B、BsB2
第２要素：BsI、BsI2、BsI3
ここで、第１要素の拡散種は、比較的低速の拡散種（低速拡散種）であり、平衡状態に達していないものとして扱われる。
一方、第２要素の拡散種は、比較的高速の拡散種（高速拡散種）であり、平衡状態に達しているものとして扱われる。
【００４６】
また、各時間領域毎に、それぞれ、計算における時間刻み間隔Δｔを設定する。本実施形態では、下記の通り、すなわち図９に示すように設定する。
第１の時間領域（ｔ≦ｔ１）では、Δｔ＝ｄｔｍｉｎ（ただし、ｄｔｍｉｎは、定数）とする。この場合、Δｔは、一定である。
第２の時間領域（ｔ１＜ｔ≦ｔ２）では、Δｔ＝ＳＱＲＴ（ｔ）・ｄｔｖａｒ（ただし、ｄｔｖａｒは、定数）とする。この場合、Δｔは、アニール経過時間ｔに応じて増大する。
第３の時間領域（ｔ２＜ｔ≦ｔ３）では、Δｔ＝ＳＱＲＴ（ｔ）・ｄｔｖａｒ（ただし、ｄｔｖａｒは、定数）とする。この場合、Δｔは、アニール経過時間ｔに応じて増大する。
第４の時間領域（ｔ＞ｔ３）では、Δｔ＝ＳＱＲＴ（ｔ３）・ｄｔｖａｒ（ただし、ｄｔｖａｒは、定数）とする。この場合、Δｔは、一定である。
【００４７】
また、アニール温度Ｔに応じて、ホウ素の拡散を予測するための時間に関する各パラメーターを設定する。
具体的には、アニール温度Ｔに応じて、各時間領域の時間帯、すなわち、前記時間ｔ１、ｔ２およびｔ３をそれぞれ設定する。これにより、第１の時間領域、第２の時間領域、第３の時間領域および第４の時間領域が、それぞれ、設定される。
【００４８】
また、各時間領域毎に、それぞれ、アニール温度Ｔに応じて、時間刻み間隔Δｔを設定する。この場合、前記ｄｔｍｉｎおよびｄｔｖａｒをそれぞれ設定する。このｄｔｍｉｎおよびｄｔｖａｒの設定と、前記時間ｔ３の設定とにより、各時間領域において、それぞれ、Δｔが設定される。
また、シリコン基板中に深さ方向１次元、または、深さ、横方向２次元のホウ素濃度分布を仮定し、これらの設定値および拡散方程式等に基づいて、ホウ素の濃度分布を時間発展させる（計算する）。そして、所定のアニール時間ｔｍａｘまで計算し、最終的なホウ素の濃度分布を得る。
この場合、第１要素に振り分けられた拡散種の濃度（濃度分布）は、下記式１で示す第１の関係式（拡散方程式）を用いて算出し、第２要素に振り分けられた拡散種の濃度は、下記式２で示す第２の関係式（定常状態による式）と、前記算出した第１要素に振り分けられた拡散種の濃度とを用いて算出する。
【００４９】
【数３】
2005057143.tif 000005
【００５０】
【数４】
2005057143.tif 000006
【００５１】
このように、平衡状態に達した拡散種の濃度を、他の平衡状態に達していない拡散種の濃度から求めるとともに、時間刻み間隔Δｔを比較的大きくすることにより、計算量が少なくなり、迅速にホウ素の濃度分布を求めることができる。
また、アニール温度Ｔに応じて、時間に関する各パラメーターを設定するので、任意のアニール温度Ｔ、特に、半導体のプロセスに使用される温度領域（例えば、０〜１２００℃程度）のすべてにおいて、計算が破綻することなく、確実に、ホウ素の濃度分布を求めることができる。
【００５２】
ここで、時間刻み間隔Δｔは、計算が破綻しない範囲において、可能な限り大きい値に設定するのが好ましい。これにより、計算量をより少なくすることができ、より迅速にホウ素の濃度分布を求めることができる。
また、例えば、半導体のプロセスに使用される温度領域（例えば、０〜１２００℃程度）のすべてに有効で最適な時間に関する各パラメーターを設定することにより、シミュレーション毎に有効なパラメーターを探し、かつ最適化する必要がなくなり、プロセス開発、設計等を容易かつ迅速に行うことができる。
【００５３】
次に、アニール温度Ｔに応じた時間に関する各パラメーターの設定の一例について説明する。
本実施形態では、アニール温度Ｔが、０℃以上、７００℃未満（０℃≦Ｔ＜７００℃）の場合、７００℃以上、８００℃未満（７００℃≦Ｔ＜８００℃）の場合、８００℃以上、１０００℃未満（８００℃≦Ｔ＜１０００℃）の場合、１０００℃以上、１２００℃以下（１０００℃≦Ｔ≦１２００℃）の場合のそれぞれにおいて、時間に関するパラメーターのセットを用意する。この場合、計算が破綻せず、かつ、計算量が最も少なくなるように、パラメーターセットを用意するのが好ましい。
ここで、０℃≦Ｔ＜７００℃、７００℃≦Ｔ＜８００℃、８００℃≦Ｔ＜１０００℃、１０００℃≦Ｔ≦１２００℃の各温度領域において、計算が破綻しないための各パラメーターｔ１、ｔ２、ｔ３、ｄｔｍｉｎおよびｄｔｖａｒの範囲は、下記表１に示す通りである。
【００５４】
【表１】
2005057143.tif 000007
【００５５】
アニール時間Ｔまでの計算量をより少なくするには、各パラメーターの値をより大きくすればよい。
従って、上記表１から判るように、０℃≦Ｔ＜７００℃、７００℃≦Ｔ＜８００℃、８００℃≦Ｔ＜１０００℃、１０００℃≦Ｔ≦１２００℃の各温度領域において、計算が破綻せず、かつ計算量が最も少ないパラメーターの組み合わせは、dt_min=10^-11、t_1=10^-7、t_3=10^-2であって、
(t_2, dt_var)=(10^-2,10^-4)、(10^-3,10^-4)、(10^-4,10^-4)、(10^-4,10^-5)
である。これらを各々、時間パラメーターセット１〜４として、図１０に示す後述する使用パラメーター記憶部４に記憶する。
【００５６】
なお、時間パラメーターセット１は、dt_min=10^-11、t_1=10^-7、t_2=10^-2、t_3=10^-2、dt_var=10^-4である。
また、時間パラメーターセット２は、dt_min=10^-11、t_1=10^-7、t_2=10^-3、t_3=10^-2、dt_var=10^-4である。
また、時間パラメーターセット３は、dt_min=10^-11、t_1=10^-7、t_2=10^-4、t_3=10^-2、dt_var=10^-4である。
また、時間パラメーターセット４は、dt_min=10^-11、t_1=10^-7、t_2=10^-4、t_3=10^-2、dt_var=10^-5である。
【００５７】
図１０は、本発明の目的物質の拡散予測方法を実行するシミュレータの実施形態の主要部を示すブロック図、図１１は、図１０に示すシミュレータの制御手段の制御動作を示すフローチャートである。
図１０に示すように、シミュレータ１は、初期条件入力部５と、シミュレーション実行部３と、使用パラメーター記憶部４と、結果記憶部６と、結果出力部７とを有している。
【００５８】
初期条件入力部５は、例えば、キーボードやマウス等で構成される。シミュレーション実行部３は、例えば、ＭＰＵ（Micro Processing Unit）や、ＣＰＵ（Central Processing Unit）等の演算装置で構成される。使用パラメーター記憶部４は、例えば、ハードディスク、メモリー等で構成される。結果記憶部６は、例えば、ハードディスク、メモリー等で構成される。結果出力部７は、例えば、ディスプレイ（表示部）等で構成される。なお、前記シミュレーション実行部３、使用パラメーター記憶部４および結果記憶部６により、制御手段２が構成されている。
使用者は、初期条件入力部５から、拡散を予測する（求めたい）アニール条件として、アニール温度（温度）Ｔと、アニール時間（時間）ｔｍａｘとを、それぞれ入力する。そして、図示しない開始スイッチをオンすると、図１１に示すプログラムがシミュレーション実行部３により実行される。
【００５９】
図１１に示すように、まず、温度Ｔが７００℃未満（Ｔ＜７００℃）であるか否かを判断し（ステップＳ１０１）、温度Ｔが７００℃未満である場合には、時間パラメーターセット１を使用パラメーター記憶部４から読み出して設定し（ステップＳ１０５）、ステップＳ１０９へ進む。
また、ステップＳ１０１において温度Ｔが７００℃未満ではないと判断した場合には、温度Ｔが７００℃以上、８００℃未満（７００℃≦Ｔ＜８００℃）であるか否かを判断し（ステップＳ１０２）、温度Ｔが７００℃以上、８００℃未満である場合には、時間パラメーターセット２を使用パラメーター記憶部４から読み出して設定し（ステップＳ１０６）、ステップＳ１０９へ進む。
【００６０】
また、ステップＳ１０２において温度Ｔが７００℃以上、８００℃未満ではないと判断した場合には、温度Ｔが８００℃以上、１０００℃未満（８００℃≦Ｔ＜１０００℃）であるか否かを判断し（ステップＳ１０３）、温度Ｔが８００℃以上、１０００℃未満である場合には、時間パラメーターセット３を使用パラメーター記憶部４から読み出して設定し（ステップＳ１０７）、ステップＳ１０９へ進む。
【００６１】
また、ステップＳ１０３において温度Ｔが８００℃以上、１０００℃未満ではないと判断した場合には、温度Ｔが１０００℃以上（Ｔ≧１０００℃）であるか否かを判断し（ステップＳ１０４）、温度Ｔが１０００℃以上である場合には、時間パラメーターセット４を使用パラメーター記憶部４から読み出して設定し（ステップＳ１０８）、ステップＳ１０９へ進む。
【００６２】
次いで、アニール経過時間（時間）ｔがｔ１以下（ｔ≦ｔ１）であるか否かを判断し（ステップＳ１０９）、時間ｔがｔ１以下である場合には、ステップＳ１１３へ進む。
また、ステップＳ１０９において時間ｔがｔ１以下ではないと判断した場合には、時間ｔがｔ１を超え、ｔ２以下（ｔ１＜ｔ≦ｔ２）であるか否かを判断し（ステップＳ１１０）、時間ｔがｔ１を超え、ｔ２以下である場合には、ステップＳ１１４へ進む。
【００６３】
また、ステップＳ１１０において時間ｔがｔ１を超え、ｔ２以下ではないと判断した場合には、時間ｔがｔ２を超え、ｔ３以下（ｔ２＜ｔ≦ｔ３）であるか否かを判断し（ステップＳ１１１）、時間ｔがｔ２を超え、ｔ３以下である場合には、ステップＳ１１５へ進む。
また、ステップＳ１１１において時間ｔがｔ２を超え、ｔ３以下ではないと判断した場合には、時間ｔがｔ３を超えている（ｔ＞ｔ３）か否かを判断し（ステップＳ１１２）、時間ｔがｔ３を超えている場合には、ステップＳ１１６へ進む。
【００６４】
ステップＳ１１３では、計算における時間刻み間隔Δｔ＝ｄｔｍｉｎ、第１要素の集合として、下記Ｉ１の拡散種、第２要素の集合として、下記Ｊ１の拡散種を、それぞれ、設定する。
Ｉ１：BsI、BsI2、BsI3、BsB、BsB2、BsB3、BsBI、BsBI2、Bs2B、BsB2
Ｊ１：なし
ステップＳ１１４では、Δｔ＝ＳＱＲＴ（ｔ）・ｄｔｖａｒ、第１要素の集合として、下記Ｉ２の拡散種、第２要素の集合として、下記Ｊ２の拡散種を、それぞれ、設定する。
Ｉ２：BsI2、BsI3、BsB、BsB2、BsB3、BsBI、BsBI2、Bs2B、BsB2
Ｊ２：BsI
ステップＳ１１５では、Δｔ＝ＳＱＲＴ（ｔ）・ｄｔｖａｒ、第１要素の集合として、下記Ｉ３の拡散種、第２要素の集合として、下記Ｊ３の拡散種を、それぞれ、設定する。
Ｉ３：BsB、BsB2、BsB3、BsBI、BsBI2、Bs2B、BsB2
Ｊ３：BsI、BsI2、BsI3
ステップＳ１１６では、Δｔ＝ＳＱＲＴ（ｔ３）・ｄｔｖａｒ、第１要素の集合として、下記Ｉ４の拡散種、第２要素の集合として、下記Ｊ４の拡散種を、それぞれ、設定する。
Ｉ４：BsB、BsB2、BsB3、BsBI、BsBI2、Bs2B、BsB2
Ｊ４：BsI、BsI2、BsI3
【００６５】
次いで、ステップＳ１１７、ステップＳ１１８を実行する。すなわち、前述したように、第１要素に振り分けられた拡散種の濃度を、上記式１で示す第１の関係式を用いて算出し、その結果を図示しないメモリーに記憶（記憶済みの場合は、更新）し、第２要素に振り分けられた拡散種の濃度を、上記式２で示す第２の関係式と、前記算出した第１要素に振り分けられた拡散種の濃度とを用いて算出、その結果を図示しないメモリーに記憶（記憶済みの場合は、更新）する。
【００６６】
次いで、時間ｔに時間刻み間隔Δｔを加算する（ｔ＝ｔ＋Δｔ）（ステップＳ１１９）。
次いで、時間ｔがｔｍａｘを超えている（ｔ＞ｔｍａｘ）か否かを判断し（ステップＳ１２０）、時間ｔがｔｍａｘを超えていない場合には、ステップＳ１０９へ戻り、再度、ステップＳ１０９以降を実行する。
【００６７】
一方、ステップＳ１２０において時間ｔがｔｍａｘを超えていると判断した場合、すなわち、アニール時間ｔｍａｘに達した場合には、結果記憶部６に、求まったホウ素の濃度分布を記憶し、結果出力部７に、その濃度分布を表示（出力）し、このプログラムを終了する。
以上説明したように、本実施形態の目的物質の拡散予測方法によれば、任意のアニール温度Ｔ、特に、半導体のプロセスに使用される温度領域（例えば、０〜１２００℃程度）のすべてにおいて、計算が破綻することなく、精度良く、迅速かつ確実に、ホウ素の濃度分布を予測することができる。
【実施例】
【００６８】
次に、本発明の具体的実施例（実験例）について説明する。
（実施例１）
前述した時間パラメーターセット２、すなわち、dt_min=10^-11、t_1=10^-7、t_2=10^-3、t_3=10^-2、dt_var=10^-4に設定し、前述した式１で示す第１の関係式および式２で示す第２の関係式を用いて、前述したように、アニール温度Ｔ＝７００℃におけるシリコン基板中の深さ方向の各クラスターそれぞれの濃度分布（密度分布）を計算した。
この結果、計算が破綻することなく、各クラスターそれぞれの濃度分布が得られた。この計算結果（後述する比較例１に対応するもの）は、下記表２に示す通りである。
【００６９】
（比較例１）
dt_min=10^-11、t_1=10^-7、t_2=10^-4、t_3=10^-2、dt_var=10^-4（アニール温度Ｔによらず、一定の値）に設定した他は、前記実施例１と同様にして、アニール温度Ｔ＝７００℃におけるシリコン基板中の深さ方向の各クラスターそれぞれの濃度分布（密度分布）を計算した。
この結果、クラスターBsB2およびBsBI2については、それぞれ、t=1.026×10^-4秒で計算が破綻した（解が負になった）。この計算が破綻したクラスターBsB2およびBsBI2についての計算結果は、下記表２に示す通りである。
【００７０】
（実施例２）
前述した時間パラメーターセット４、すなわち、dt_min=10^-11、t_1=10^-7、t_2=10^-4、t_3=10^-2、dt_var=10^-5に設定し、前述した式１で示す第１の関係式および式２で示す第２の関係式を用いて、前述したように、アニール温度Ｔ＝１１００℃におけるシリコン基板中の深さ方向の各クラスターそれぞれの濃度分布（密度分布）を計算した。
この結果、計算が破綻することなく、各クラスターそれぞれの濃度分布が得られた。この計算結果（後述する比較例２に対応するもの）は、下記表２に示す通りである。
【００７１】
（比較例２）
dt_min=10^-11、t_1=10^-7、t_2=10^-4、t_3=10^-2、dt_var=10^-4（アニール温度Ｔによらず、一定の値）に設定した他は、前記実施例２と同様にして、アニール温度Ｔ＝１１００℃におけるシリコン基板中の深さ方向の各クラスターそれぞれの濃度分布（密度分布）を計算した。
この結果、クラスターBsIおよびBsB2については、それぞれ、t=2.502×10^-4秒で計算が破綻した（解が負になった）。この計算が破綻したクラスターBsIおよびBsB2についての計算結果は、下記表２に示す通りである。
【００７２】
【表２】
2005057143.tif 000008
【００７３】
これらの結果から、本発明によれば、任意のアニール温度Ｔにおいて、計算が破綻することなく、精度良く、迅速かつ確実に、ホウ素の濃度分布を予測することができることが判る。
これに対し、比較例では、いずれのアニール温度Ｔであっても一定（同一）の時間パラメーターを設定するので、アニール温度Ｔによっては、途中で計算が破綻してしまい、ホウ素の濃度分布を予測することができない。
【００７４】
以上、本発明の目的物質の拡散予測方法を、図示の実施形態に基づいて説明したが、本発明はこれに限定されるものではなく、各部の構成は、同様の機能を有する任意の構成のものに置換することができる。また、本発明に、他の任意の構成が付加されていてもよい。
また、前記実施形態では、結晶基板として単結晶シリコン基板、目的物質としてホウ素を、それぞれ、例に挙げて説明したが、本発明では、これらに限らず、例えば、目的物質としては、リン、ヒ素等が挙げられる。
【００７５】
また、本発明は、例えば、デジタル電子回路またはコンピュータハードウェア、ファームウェア、ソフトウェア等、またはこれらの組み合わせにおいて実現できる（適用することができる）。また、本発明は、例えば、プログラマブルプロセッサにより実行できる機械判読可能な記憶装置において具現化されるコンピュータプログラム製品においても実現できる。また、本発明は、例えば、入力データで動作し、出力を生成することによって、本発明の機能を実現する命令プログラムを実行するプログラマブルプロセッサによっても行うことができる。本発明は、データ記憶系、少なくとも一つの入力装置、少なくとも一つの出力装置からデータと命令を受け、またこれらにデータと命令を送信するよう接続された少なくとも一つのプログラマブルプロセッサを備えるプログラマブル系上で実行される一つまたは複数のコンピュータプログラムにおいて有利に実現される。各コンピュータプログラムは、例えば、高レベル手順またはオブジェクト指向のプログラム言語あるいは、希望があればアセンブリまたは機械言語において実現でき、いずれの場合も、言語はコンパイルまたは解釈された言語とすることができる。適当なプロセッサの例としては、汎用および特定用途マイクロプロセッサがある。一般に、プロセッサは、リードオンリーメモリまたはランダムアクセスメモリから命令とデータを受信する。一般に、コンピュータは、データファイルを記憶する一つまたは複数の大量記憶装置を備え、これらの装置は内蔵ハードディスクと脱着式ディスク等の磁気ディスク、磁気光ディスク、光ディスク等を備える。コンピュータプログラムの命令とデータを具現化するのに適した記憶装置には、あらゆる形態の不揮発性メモリが含まれ、その一例を挙げれば、ＥＰＲＯＭ、ＥＥＰＲＯＭ、フラッシュメモリ等の半導体メモリ装置、内蔵ハードディスクと脱着式ディスク等の磁気ディスク、光磁気ディスク、ＣＤ−ＲＯＭディスク、ＤＶＤディスク等がある。上記はいずれも、特に設計されたＡＳＩＣ（特定用途集積回路）で補足してもよく、またはその中に組み込んでもよい。
【図面の簡単な説明】
【００７６】
【図１】計算機シミュレーションを用いたプロセス設計の一般的アプローチを示す図である。
【図２】電子装置プロセスの一般的なマルチスケールのモデル化方法を示すフローチャートである。
【図３】超極浅接合プロセスのマルチスケールのモデル化方法の適用例を示す図である。
【図４】シリコン中のホウ素ホウ素対拡散のための予測反応経路に沿った構造を示す図である。
【図５】（ａ）Ｂｉ拡散、（ｂ）Ｂｓ拡散による最も安定したホウ素クラスター、Ｂｓ−Ｂｓ＝Ｂｉに関する予測メカニズムを示す図である。
【図６】Ｂ−Ｓｉｉクラスターの予測クラスター化および分解を示す図であり、（ａ）Ｂｓ−Ｓｉｉ、（ｂ）Ｂｓ−Ｓｉｉ，２、（ｃ）Ｂｓ−Ｓｉｉ，３、（ｄ）Ｂｓ−Ｓｉｉ，４、（ｅ）Ｓｉｉ，２、（ｆ）Ｓｉｉ，３である。
【図７】キネティックモンテカルロシミュレーションの略図である。
【図８】ＣＢ＝ＣＳｉｉ＝１０１９ｃｍ−３、Ｔ＝９００Ｋのキネティックモンテカルロシミュレーションの結果を示す図であり、（ａ）ｔ＝０秒、（ｂ）ｔ＝１０−６秒、（ｃ）ｔ＝１０−３秒、（ｄ）ｔ＝１００秒である。
【図９】各時間領域と計算における時間刻み間隔Δｔとの関係を示すグラフである。
【図１０】本発明の目的物質の拡散予測方法を実行するシミュレータの実施形態の主要部を示すブロック図である。
【図１１】図１０に示すシミュレータの制御手段の制御動作を示すフローチャートである。
【符号の説明】
【００７７】
１……シミュレータ  ２……制御手段  ３……シミュレーション実行部  ４……使用パラメーター記憶部  ５……初期条件入力部  ６……結果記憶部  ７……結果出力部  Ｓ１０１〜Ｓ１２０……ステップ  ２１０〜２４０……ステップ  ３１０〜３３０……ステップ  ７１０〜７３０……ステップ
【図１】
2005057143.tif 000011
【図２】
2005057143.tif 000012
【図３】
2005057143.tif 000013
【図４】
2005057143.tif 000014
【図５】
2005057143.tif 000015
【図６】
2005057143.tif 000016
【図７】
2005057143.tif 000017
【図８】
2005057143.tif 000018
【図９】
2005057143.tif 000019
【図１０】
2005057143.tif 000020
【図１１】
2005057143.tif 000021