イオン注入分布発生方法及びシミュレータ
【課題】 イオン注入分布発生方法及びシミュレータに関し、近似を効果的に用いることによって計算時間を大幅に短縮する。
【解決手段】 半導体に注入するイオンの飛程の射影Ｒp の２次の項まで考慮した射影〈Ｒp （Ｅ）〉(2) を、１次の項まで考慮した既知の射影を〈Ｒp （Ｅ）〉(1) とした時に、摂動項Δp (2) （Ｅ）を用いて、
〈Ｒp （Ｅ）〉(2) ＝〈Ｒp （Ｅ）〉(1) ＋Δp (2) （Ｅ）
とした近似式を用いた２次の摂動モデルを用いて求める。
【解決手段】 半導体に注入するイオンの飛程の射影Ｒp の２次の項まで考慮した射影〈Ｒp （Ｅ）〉(2) を、１次の項まで考慮した既知の射影を〈Ｒp （Ｅ）〉(1) とした時に、摂動項Δp (2) （Ｅ）を用いて、
〈Ｒp （Ｅ）〉(2) ＝〈Ｒp （Ｅ）〉(1) ＋Δp (2) （Ｅ）
とした近似式を用いた２次の摂動モデルを用いて求める。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明はイオン注入分布発生方法及びシミュレータに関し、Ｍｏｎｔｅ Ｃａｒｌｏ ｓｉｍｕｌａｔｏｒとほぼ同じ精度で、より短時間でイオン注入分布を発生させるための手法に関するものである。
【背景技術】
【０００２】
シリコン集積回路装置において、シリコン基板への不純物の導入はイオン注入で行われるのが一般的である。
このようなシリコン集積回路装置のプロセス構築に際しては、必要な素子構造を得るためのイオン注入条件を決定する必要があるが、このようなイオン注入条件をシミュレーションにより決定することが行われている。
非晶質層へのイオン注入分布を理論的に予想する手段としてＭｏｎｔｅ Ｃａｒｌｏ がある。これは、入射イオンと基板との相互作用を、核阻止能及び電子阻止能の物理に基づいて、入射イオンの軌跡を追跡していくものである。
【０００３】
この理論は、任意のイオンを任意の基板にイオン注入した場合の一般的な場合にも有効であり、電子阻止能をチューニングすればその精度をさらに向上させることができる。
図１５は、計算モデルであり、質量数Ｍ1 ，原子番号Ｚ1 ，エネルギーＴ1i（速度ｖ1i）のイオンが、基板を構成する質量数Ｍ2 ，原子番号Ｚ2 の原子と相互作用して伝達するエネルギーＴ2fは、散乱角度をΦ、相互作用後のイオンの速度をν1i、基板原子の速度をν2iとすると、
Ｔ2f／Ｔ1i＝２Ｍ2 ν2i2 ｓｉｎ2 （Φ／２）／〔（１／２）Ｍ1 ν1i2 〕
＝（４Ｍ2 ／Ｍ1 ）｛〔Ｍ1 ｖ1i／（Ｍ1 ＋Ｍ2 ）〕2 ／ν1i2 ｝
×ｓｉｎ2 （Φ／２）
＝〔４Ｍ2 Ｍ1 ／（Ｍ1 ＋Ｍ2 ）2 〕ｓｉｎ2 （Φ／２）・・・（１）
と表現される。
【０００４】
ここで、イオンと基板原子の距離をｒ、衝突パラメータをｂ、ポテンシャルエネルギーをＶ（ｒ）とすると、伝達エネルギーＴ2fは、下記の式（２）として求まる。
【数１】
つまり、注入されたイオンは、核との相互作用により、
ΔＥn ＝Ｔ2f
のエネルギーを失う。
【０００５】
また、相互作用に伴う散乱角度Φは、下記の式（３）で表される。
【数２】
ここで、半径ｂの円周上の位置、即ち角度θは、Ｒａｎｄ（ｎ）をｎが０から１の間の乱数とすると、
θ＝２πＲａｎｄ（ｎ）
の関係から求まる。
【０００６】
これにより、衝突後のエネルギーと方向が決まる。
次に、新たに注入エネルギーで同様の計算を繰り返し、イオンの軌跡をトレースしていけば良い。
しかし、実際の衝突のたびに、上記の式（３）で表される積分を毎回実行するのは計算コストが膨大になるため、ＺｉｅｇｌｅｒはＭａｇｉｃ ｆｏｒｍｕｌａを提案し、この積分を簡単なプロトコルで解いている。
また、シミュレータＴＳＵＰＲＭでは、計算結果をテーブル化しその内挿でこの値を評価している。
【０００７】
なお、上記の式（３）における各パラメータを規格化したユニバーサルな変数で表すと、下記の式（４）で表される。
【数３】
但し、ａB をボーア径とすると、
ρ＝ｒ／ａU ，
η＝ｂ／ａU ，
ａU ＝０．８８５４ａB ／（Ｚ10.23 ＋Ｚ20.23 ）
である。
また、エネルギーは、
ε＝ＥC ／（Ｚ1 Ｚ2 ｅ2 ／ａU ）
Ｅc ＝（１／２）Ｍc ｖ1i2
１／Ｍc ＝１／Ｍ1 ＋１／Ｍ2
である。
また、ポテンシャルエネルギーＶ（ｒ）として、下記のＺｉｅｇｌｅｒ−Ｌｉｔｍａｒｋ−Ｂｉｅｒｓａｋ（ＺＬＢ）のポテンシャルエネルギーを用いる（例えば、非特許文献１参照）。
Ｖ（ｒ）＝（ｅ2 Ｚ1 Ｚ2 ／ｒ）ｆ（ρ）
但し、
【数４】
【０００８】
しかしながら、上述の計算コスト低減の手法を導入しても、Ｍｏｎｔｅ Ｃａｒｌｏシミュレーションは、粒子の各軌跡を追うため、統計誤差を減らすためには数万個以上の計算をする必要があり、時間がかかるという問題がある。
【０００９】
そこで、電子阻止能Ｓe 、核阻止能Ｓn が与えられた時、イオン注入分布が従うべき積分方程式が、Ｌｉｎｄｈａｒｔ，Ｓｃｈａｒｆ，Ｓｃｈｉｏｔｔによって提案され（例えば、非特許文献２参照）、このモデルはＬＳＳ理論と呼ばれている。
【００１０】
この理論では、粒子の軌跡を追跡することなしに注入条件が決まれば、積分方程式を解くことによって分布の任意の次数の分布モーメントのエネルギー依存性までを即座に計算できる。
【００１１】
この場合、積分方程式を展開し、微分方程式に還元し、それを解いていく必要があり、そのステップで近似が入る。
即ち、伝達エネルギーが初期エネルギーに比べて小さいと仮定してモーメントをエネルギーに関してＴａｙｌｏｒ展開するが、この展開する次数が解くべき微分方程式の階数となる。
【００１２】
また、散乱角度もエネルギーに関してＴａｙｌｏｒ展開するが、これは解くべき方程式の係数と関連し階数とは無関係である。解析モデルは一階の線形微分方程式の場合のみ任意の係数で一般的に得られる。
このＬＳＳ理論ではモーメントに関しては１次、散乱角度に関しては２次までＴａｙｌｏｒ展開し、線形１階の微分方程式を解きＲp 及びΔＲp の解析モデル式を導出している。
【００１３】
ここで、図１６を参照して、ＬＳＳ理論における注入イオンの飛程Ｒに関する積分方程式を説明する。
図１６はイオンの飛程Ｒの模式図であり、エネルギーＥで注入されたイオンは基板原子と相互作用しながらエネルギーを失い方向を変えながら図に示すように進んでいき、全エネルギーを失い基板中で静止する。
【００１４】
この時、イオンが進んだ軌跡の線積分を飛程Ｒとし、飛程Ｒの注入方向への射影をＲp 、横方向への広がりをＲT とする。
また、横方向への広がりのｘ成分をΔｘ、ｙ成分をΔｙとする。
なお、横方向への広がりは、図面及び以降の積分方程式或いは微分方程式においてはＲに「垂直記号」のサフィックスを付けた記号で表すが、明細書本文中では、明細書作成の都合上、「ＲT 」で表す。
【００１５】
表面に垂直にエネルギーＥでイオン注入されたイオンがＲとＲ＋ｄＲの間で止まる確率をＰ（Ｅ，Ｒ）とする。この場合、
【数５】
である。
また、Ｒに関してｍ次のモーメントを
【数６】
と定義する。
【００１６】
イオンが表面からｄＲ進む間に原子および電子と相互作用する確率は、
ＮｄＲ∫ｄσn ＋ＮｄＲ∫ｄσe ・・・（８）
であり、そのそれぞれの相互作用でイオンのエネルギーは、
Ｅ−Ｔn ，Ｅ−Ｔe
に減少すると仮定する。
ここで、ｄσn ，ｄσe は核阻止能及び電子阻止能と関連する微分断面積である。
【００１７】
よって、このイオンが原子および電子と相互作用し、エネルギーを失って飛程Ｒに止まる確率は、
【数７】
で与えられる。
また、ｄＲ進む間に衝突しない確率は、
１−（ＮｄＲ∫ｄσn ＋ＮｄＲ∫ｄσe ） ・・・（１０）
である。この間にエネルギーは失われないから、イオンがＲに止まる確率は、
［１−（ＮｄＲ∫ｄσn ＋ＮｄＲ∫ｄσe ）］Ｐ（Ｅ，Ｒ−ｄＲ） ・・・（１１）
となる。よって、
【数８】
となる。
【００１８】
ここで、ｄＲを無限小に持っていく極限を考えると、
【数９】
となる。両辺にＲm をかけてＲについて０から∞まで積分すると、
【数１０】
となる。
【００１９】
式（１４）の左辺を部分積分し、また右辺の積分の順序を入れ換えると、
【数１１】
式（１５）における左辺の第１項は０である。また
【数１２】
と表現するが、これはＲm の期待値に相当する。
【００２０】
ここで、式（１６）を用いると、式（１５）は、
【数１３】
となり、これが飛程Ｒに関するｍ次の積分方程式である。
ここで、飛程Ｒに関しては１次に関してのみ解析する。式（１７）においてｍ＝１とおいて
【数１４】
を得る。これが、飛程Ｒの従うべき積分方程式である。
【００２１】
次に、式１７をＴn ，Ｔe ≪Ｅの近似のもとで１次の項までを残し解いていく。この近似は伝達エネルギーが小さいということを仮定をしている。
この場合、核との相互作用に関しては広角散乱を無視する近似に相当し、電子との相互作用に関しては、エネルギーが高いほどＴe は大きくなる。このため、比較的エネルギーの低い場合に近似の精度は良くなってくる。
この解を１次のものであることを意識して〈Ｒ（Ｅ）〉(1) と表現する。
【００２２】
次に、式（１８）の両辺をＮで割るとともにＴａｙｌｏｒ展開することにより、
【数１５】
となる。ここで、
Ｓn ＝∫Ｔn ｄσn ，Ｓe ＝∫Ｔe ｄσe ・・・（２０）
を利用している。式（１９）より、良く知られている、
【数１６】
が導出される。
【００２３】
次に、図１７を参照して、Ｒp （Ｅ，ｃｏｓφ），ＲT （Ｅ，ｃｏｓφ），Ｒp （Ｅ），ＲT （Ｅ）の幾何学的関係を説明する。
図１７は、Ｒp （Ｅ，ｃｏｓφ），ＲT （Ｅ，ｃｏｓφ），Ｒp （Ｅ），ＲT （Ｅ）の幾何学的関係の説明図であり、エネルギーＥ、角度φで入射したイオンが基板内のＡ点に静止した状況を模式図的に示している。
入射角度０の軸に垂直な面上での横方向の広がりをＲT （Ｅ，ｃｏｓφ）とすると、入射方向に垂直な軸に対する射影Ｒp （Ｅ）が、入射方向に垂直な面に対する横方向広がりがＲT （Ｅ）と考えることができる。よって、
【数１７】
【００２４】
以上の幾何学的関係を基にして、次に、Ｒp ，ＲT ，Ｒp ＲT の従うべき積分方程式を導出していく。
エネルギーＥで、角φで散乱されたイオンの射影飛程がＲp とＲp ＋ｄＲp の間で止まる確率をＰp （Ｅ，Ｒp ，ｃｏｓφ）とする。
φ＝０の場合の確率をＰp （Ｅ，Ｒp ）とする。また、Ｒp に関してｍ次のモーメントを、
【数１８】
と定義する。
【００２５】
入射イオンがｄＲp 進む間に原子および電子と相互作用する確率は、飛程Ｒの場合と同様に、
ＮｄＲp ∫ｄσn ＋ＮｄＲp ∫ｄσe ・・・（２５）
であり、この相互作用でのエネルギーは、
Ｅ−Ｔn ，Ｅ−Ｔe
に減少する。
【００２６】
この場合、電子との相互作用による散乱角は０とみなすので、このイオンが原子および電子と相互作用し、エネルギーＴn ，Ｔe を失ってＲp に止まる確率は、
【数１９】
で与えられる。
また、ｄＲ進む間に衝突しない確率は、
１−（ＮｄＲp ∫ｄσn ＋ＮｄＲp ∫ｄσe ） ・・・（２７）
である。この間にエネルギーは失われないから、イオンがＲp に止まる確率は、
［１−（ＮｄＲp ∫ｄσn ＋ＮｄＲp ∫ｄσe ）］Ｐ（Ｅ，Ｒp −ｄＲp ）
・・・（２８）
となる。よって、
【数２０】
となる。
【００２７】
ここで、ｄＲp を無限小にもっていく極限を考えると、
【数２１】
となり、両辺にＲm p をかけてＲp について０から∞まで積分すると、
【数２２】
とｍ次の〈Ｒm p （Ｅ）〉に関する積分方程式が導出される。
【００２８】
次に、分布の横方向への広がりを表すＲT について考える。
エネルギーＥ、入射角φで注入されたイオンがＲT とＲT ＋ｄＲT の間に止まる確率をＰT （Ｅ，ＲT ，ｃｏｓφ）とする。
この時、入射角φのイオンがｄＲT 進む間に原子および電子と相互作用する確率は、
【数２３】
である。よって、このイオンが相互作用しエネルギーＴn ，Ｔe を失ってＲT に止まる確率は、
【数２４】
で与えられる。
【００２９】
この場合も電子との相互作用では角度の変化はなく、さらにこの場合はＲT も変化しないことを仮定している。また、ｄＲT 進む間に衝突しない確率は、
【数２５】
である。この間にエネルギーは失われず、角度も変化しないためＲT も変化しない。よって、イオンがＲT に止まる確率は、
【数２６】
となる。
【００３０】
以上より、
【数２７】
となる。
【００３１】
ここで、ｄＲT を無限小にもっていく極限を考えると、
【数２８】
となり、両辺にＲm T をかけてＲT について０から∞まで積分すると、
【数２９】
と〈Ｒm T （Ｅ）〉に関する積分方程式が導出される。
【００３２】
まず、飛程の射影Ｒp を求めるために、式（３１）においてｍ＝１とすると下記の式（３９）となる。
【数３０】
ここで、１次のテーラー展開した解を〈Ｒp （Ｅ）〉(1) とおくと、下記の式（４０）となる。
【数３１】
【００３３】
これから、下記の式（４１）が得られる。
【数３２】
但し、α（Ｅ）及びα（Ｅ′）は、下記の式（４２）で表される。
【数３３】
【００３４】
次に、２次のモーメントΔＲp 及びΔＲptを求めるために、式（３１）及び式（３８）においてｍ＝２とすると下記の式（４３）及び式（４４）となる。
【数３４】
【００３５】
ここで、Ｌｉｎｄｈａｒｔは下記の式（４５）及び式（４６）で定義される変数〈Ｒ2 c （Ｅ）〉及び〈Ｒ2 r （Ｅ）〉を導入した。
【数３５】
【００３６】
〈ΔＲ2 p （Ｅ）〉と〈Ｒ2 T （Ｅ）〉は、変数〈Ｒ2 c （Ｅ）〉及び〈Ｒ2 r （Ｅ）〉と下記の式（４７）及び式（４８）で関連づけられる。
【数３６】
【００３７】
また、横方向の偏差〈ΔＲ2 pt（Ｅ）〉は〈ΔＸ2 〉と関連し、〈Ｒ2 T （Ｅ）〉そのものではない。したがって、
〈Ｒ2 T （Ｅ）〉＝〈ΔＸ2 〉＋〈ΔＸ2 〉＝２〈ΔＲ2 pt（Ｅ）〉・・・（４９）
となり、これより、
〈ΔＲ2 pt（Ｅ）〉＝〈Ｒ2 T （Ｅ）〉／２ ・・・（５０）
となる。
【００３８】
Ｌｉｎｄｈａｒｔはこれらの変数を導入することにより、下記の式（５１）及び式（５２）の独立な微分方程式を得ることに成功した。
【数３７】
【００３９】
また、Ｇｉｂｂｏｎｓは、このモデル式にエネルギーに依存しない核阻止能を代入し、Ｓｉ中のＢ，Ｎ，Ａｌ，Ｐ，Ｇａ，Ａｓ，Ｉｎ，ＳｂのＲp ，ΔＲp を評価し、いくつかの実験データと比較し、良い一致を得ている。
さらに、３次までのモーメントを数値的に評価することが検討されている。
【非特許文献１】Ｊ．Ｆ．Ｚｉｅｇｌｅｒ，Ｊ．Ｐ．Ｂｉｅｒｓａｃｋ，ａｎｄ Ｕ．Ｌｉｔｔｍａｒｋ，Ｔｈｅ ｓｔｏｐｐｉｎｇ ａｎｄ ｒａｎｇｅ ｏｆ ｉｏｎｓ ｉｎ ｓｏｌｉｄ，Ｐｅｒｇａｍｏｎ，１８８５
【非特許文献２】Ｊ．Ｌｉｎｄｈａｒｔ，Ｍ．Ｓｃｈａｒｆｆ，Ｈ．Ｓｃｈｉｏｔｔ，Ｍａｔ．Ｆｔｓ．Ｍｅｄｄ．Ｖｉｄ．Ｓｃｌｓｋ，ｖｏｌ．３３，ｐｐ．１−３９，１９６３
【発明の開示】
【発明が解決しようとする課題】
【００４０】
しかし、ＬＳＳ理論では、積分方程式を解く際に近似が必要で、これまでに提案されているモデルではＲp は問題ないが、ΔＲp は格段に精度が落ちるか、計算できないという問題がある。
【００４１】
また、上述のように、３次までのモーメントを数値的に評価することが検討されているが、これらの数値的に解析したデータは系統的に比較されていない。また、比較するうえでも、３次までのモーメントだけでは現在利用されているＰｅａｒｓｏｎ分布を利用することが出来ないという問題点もある。
【００４２】
即ち、核阻止能モデル、電子阻止能モデルの種類、積分方程式の展開の仕方によってＬＳＳ理論から導かれる結果はさまざまになる。ＬＳＳ理論で扱われている近似の精度は同じ核阻止能モデル、電子阻止能モデルを使ったＭｏｎｔｅ Ｃａｒｌｏ計算と比較することによって正確に評価される。
【００４３】
しかしながら、精度の高い核阻止能、電子阻止能を使ってＭｏｎｔｅ ＣａｒｌｏとＬＳＳ理論の比較がこれまでなされていない。さらに実際の分布は４次のモーメントまで利用して初めて精度よく表現される。
したがって、これまでのＬＳＳ理論では、現在利用されているＰｅａｒｓｏｎ関数と対応づけることはできない。
【００４４】
また、衝突パラメータｂと散乱角度Φの相関からイオン注入分布を求める場合も、Ｍｏｎｔｅ Ｃａｒｌｏは粒子の各軌跡を追うため、統計誤差を減らすためには数万個以上の計算をする必要があるため時間がかかる。
この計算を進めていく上で一番時間のかかるのが散乱角度を評価するための積分である。
【００４５】
また、衝突と次の衝突の間の距離Ｌをいくらにとるのかも計算時間に影響をおよぼす。
これまでは平均原子間距離Ｌmin が利用されており、この場合は相互作用が無視できるような場合でも散乱を計算してしまうため計算時間が膨大になるという問題がある。
【００４６】
したがって、本発明は、実験データを再現するＭｏｎｔｅ Ｃａｒｌｏ ｓｉｍｕｌａｔｏｒと同じ核阻止能、電子阻止能を用いたＬＳＳ理論の積分方程式から、２次、３次の摂動計算によってＲp ，ΔＲp ，ΔＲpt，γの解析モデルを導出して計算時間を大幅に短縮することを目的とする。
【００４７】
また、Ｍｏｎｔｅ Ｃａｒｌｏを用いる場合にも、散乱角度を評価するための積分を解析的に表現する近似式で置き換えて計算時間を大幅に短縮することを目的とする。
【課題を解決するための手段】
【００４８】
本発明の一観点によると、半導体に注入するイオンの飛程の射影Ｒp の２次の項まで考慮した射影〈Ｒp （Ｅ）〉(2) を、１次の項まで考慮した既知の射影を〈Ｒp （Ｅ）〉(1) とした時に、摂動項Δp (2) （Ｅ）を用いて、
〈Ｒp （Ｅ）〉(2) ＝〈Ｒp （Ｅ）〉(1) ＋Δp (2) （Ｅ）
とした近似式を用いた２次の摂動モデルを用いて求めるイオン注入発生方法が提供される。
【００４９】
また、別の観点によると、基板に注入したイオンの散乱角度Φをモンテ・カルロ法で求める際に、散乱角度Φを、飽和する一定高角度の項、低角度の項、中間角度の項の３つの項の和の逆数で解析式として表すイオン注入分布発生方法が提供される。
【００５０】
さらに、別の観点によると、上述の各種のイオン注入分布発生方法に関する計算式が組み込まれたシミュレータが提供される。
【発明の効果】
【００５１】
本発明によれば、精度の高い解析モデルを導出し、さらに高次のモーメントを高精度で予想できるようになった。これにより、現在標準的に用いられているＰｅａｒｓｏｎ分布のパラメータを初めて解析的に求めることが可能になる。
したがって、任意不純物を任意非晶質基板にイオン注入した分布のエネルギー依存性が瞬時に予想できるようになる。
また、この分布発生ソフトにより横方向の分布の深さ依存性も予想できる。
【００５２】
さらに、Ｍｏｎｔｅ Ｃａｒｌｏを用いる場合にも、散乱角度を評価するための積分を解析的に表現する近似式で置き換えてることによって、計算時間を大幅に短縮することが可能になる。
【発明を実施するための最良の形態】
【００５３】
まず、１次乃至３次の各阻止能について、計算しておく。
１次の核阻止能Ｓn （Ｅ）は、下記の式（５３）で与えられる。
【数３８】
また、ユニバーサルな核阻止能Ｓn （ε）は、下記の式（５４）で与えられる。
【数３９】
【００５４】
ここで、上述のＺＬＢポテンシャルを利用すると、ユニバーサルな核阻止能Ｓn （ε）は、下記の式（５５）で近似される。
【数４０】
【００５５】
この阻止能は高次のものにも拡張され、２次のモーメントに対応する拡張された核阻止能Ωn 2 は、下記の式（５６）で表される。
【数４１】
但し、ユニバーサルな拡張された核阻止能Ωn 2 （ε）は下記の式（５７）で表される。
【数４２】
これも、下記の式（５８）で近似的に表現される。
【数４３】
【００５６】
さらに、３次のモーメントに対応する拡張された核阻止能Λn 3 は、下記の式（５９）で表される。
【数４４】
但し、ユニバーサルな拡張された核阻止能Λn 3 （ε）は下記の式（６０）で表される。
【数４５】
これも、下記の式（６１）で近似的に表現される。
【数４６】
【００５７】
ここで、上記式（２２）及び式（２３）に基づいて、以下の計算で利用する各パラメータＲp 2 （Ｅ，ｃｏｓφ），Ｒp 3 （Ｅ，ｃｏｓφ）、Ｒp ＲT 2 （Ｅ，ｃｏｓφ）を予め下記の通り計算しておく。
【数４７】
【００５８】
ここで、上記の図１７に示したように、空間の対称性から点Ｃの存在する平面上でＣを中心とする半径ＲT （Ｅ）の円上の点はＡと等価と考えることができる。
よって、βは０から２πの値を同じ確率で取ると仮定できる。
そこで、式（６２）乃至式（６４）をβに関して０から２πまで積分したものを２πで割りその平均値を評価する。
【００５９】
【数４８】
を参考にしてこれらの平均値は、
【数４９】
となる。ここで、
【数５０】
である。
【００６０】
次に、γ等の高次のモーメントを求めるために必要なクロスタームに関する積分方程式を下記の式（７３）として導出した。
【数５１】
が得られる。
【００６１】
次に、飛程の射影Ｒp について、上記の式（３９）を２次まで展開した解を〈Ｒp （Ｅ）〉(2) とおくと、
【数５２】
となる。
但し、Ωn 2 ＝∫Ｔn 2 ｄσn ，Ωe 2 ＝∫Ｔe 2 ｄσe
を利用している。
なお、一般に、Ωe 2 ≪Ωn 2 であるので、実際の計算ではΩe 2 は無視する。
【００６２】
上記の式（７４）は２階の微分方程式になり、解析的に解くことはできない。そこで摂動近似で２次の近似式〈Ｒp （Ｅ）〉(2) を、計算可能な既知の〈Ｒp （Ｅ）〉(1) と摂動項Δp (2) を用いて、
〈Ｒp （Ｅ）〉(2) ＝〈Ｒp （Ｅ）〉(1) ＋Δp (2) （Ｅ） ・・・（７５）
と表現する。
【００６３】
この式（７５）を上記の式（７４）に代入し、Ω2 とΔp (2) の積は３次の微小項と見なして落としたのち、上記の１次展開に関する式（４０）を利用することによって、
【数５３】
が得られ、これより、摂動項Δp (2) は、下記の式（７７）で表される。
【数５４】
但し、式（７７）におけるζp (2) （Ｅ）は下記の式（７８）で表される。
【数５５】
【００６４】
このように、未知の〈Ｒp （Ｅ）〉(2) を、計算可能な既知の〈Ｒp （Ｅ）〉(1) と摂動項Δp (2) により近似することによって、飛程の射影Ｒp の２次の摂動モデルによる値〈Ｒp （Ｅ）〉(2) を数値的に計算することが可能になる。
【００６５】
同様に、飛程の射影Ｒp について、上記の式（３９）を３次まで展開した解を〈Ｒp （Ｅ）〉(3) とおくと、
【数５６】
となる。
但し、Λn 3 ＝∫Ｔn 3 ｄσn ，Λe 3 ＝∫Ｔe 3 ｄσe
を利用している。
【００６６】
次に、〈Ｒp （Ｅ）〉(3) はそのままでは計算できないので、上記の式（７７）により既知になった〈Ｒp （Ｅ）〉(2) と摂動項Δp (3) を用いて、
〈Ｒp （Ｅ）〉(3) ＝〈Ｒp （Ｅ）〉(2) ＋Δp (3) （Ｅ） ・・・（８０）
と表現する。
【００６７】
この近似式を用いて２次の場合と同様に計算すると、摂動項Δp (3) は、下記の式（８１）で表される。
【数５７】
但し、式（８１）におけるζp (3) （Ｅ）は下記の式（８２）で表される。
【数５８】
【００６８】
この場合も未知の〈Ｒp （Ｅ）〉(3) を、既知となって計算可能な〈Ｒp （Ｅ）〉(2) と摂動項Δp (3) により近似することによって、飛程の射影Ｒp の３次の摂動モデルによる値〈Ｒp （Ｅ）〉(3) を数値的に計算することが可能になる。
【００６９】
次に、ΔＲp とΔＲptについて検討する。
まず、式（５１）の積分方程式を２次まで展開した解を〈Ｒc 2 （Ｅ）〉(2) とおくと、下記の式（８３）が得られる。
【数５９】
ここで、未知の〈Ｒc 2 （Ｅ）〉(2) を、計算可能な既知の〈Ｒc 2 （Ｅ）〉(1) と摂動項Δc 2(2)を用いて、
〈Ｒc 2 （Ｅ）〉(2) ＝〈Ｒc 2 （Ｅ）〉(1) ＋Δc 2(2) ・・・（８４）
とおいて、上記の式（８３）を整理すると、左辺は、
２〈Ｒp （Ｅ）〉(2) ／Ｎ＝２｛〈Ｒp （Ｅ）〉(1) ＋Δp (2) （Ｅ）｝／Ｎ であるので、下記の式（８５）となる。
【数６０】
【００７０】
したがって、摂動項Δc 2(2)は、下記の式（８６）となる。
【数６１】
但し、式（８６）におけるζc (2) （Ｅ）は下記の式（８７）で表される。
【数６２】
【００７１】
同様に、式（５１）の積分方程式を３次まで展開した解を〈Ｒc 2 （Ｅ）〉(3) とおくと、下記の式（８８）が得られる。
【数６３】
ここで、未知の〈Ｒc 2 （Ｅ）〉(3) を、上記の結果により既知となっった〈Ｒc 2 （Ｅ）〉(2) と摂動項Δc 2(3)を用いて、
〈Ｒc 2 （Ｅ）〉(3) ＝〈Ｒc 2 （Ｅ）〉(2) ＋Δc 2(3) ・・・（８９）
とおいて、上記の式（８８）を整理すると、左辺は、
２〈Ｒp （Ｅ）〉(3) ／Ｎ＝２｛〈Ｒp （Ｅ）〉(2) ＋Δp (3) （Ｅ）｝／Ｎ であるので、下記の式（９０）となる。
【数６４】
【００７２】
したがって、摂動項Δc 2(3)は、下記の式（９１）となる。
【数６５】
但し、式（９１）におけるζc (3) （Ｅ）は下記の式（９２）で表される。
【数６６】
【００７３】
次に、式（５２）の積分方程式を２次まで展開した解を〈Ｒr 2 （Ｅ）〉(2) とおくと、下記の式（９３）が得られる。
【数６７】
ここでも、未知の〈Ｒr 2 （Ｅ）〉(2) を、計算可能な既知の〈Ｒr 2 （Ｅ）〉(1) と摂動項Δr 2(2)を用いて、
〈Ｒr 2 （Ｅ）〉(2) ＝〈Ｒr 2 （Ｅ）〉(1) ＋Δr 2(2) ・・・（９４）
とおいて、上記の式（９３）を整理すると、左辺は、
２〈Ｒp （Ｅ）〉(2) ／Ｎ＝２｛〈Ｒp （Ｅ）〉(1) ＋Δp (2) （Ｅ）｝／Ｎ であるので、下記の式（９５）となる。
【数６８】
【００７４】
したがって、摂動項Δr 2(2)は、下記の式（９６）となる。
【数６９】
但し、式（９６）におけるζr (2) （Ｅ）は下記の式（９７）で表される。
【数７０】
【００７５】
同様に、式（５２）の積分方程式を３次まで展開した解を〈Ｒr 2 （Ｅ）〉(3) とおくと、下記の式（９８）が得られる。
【数７１】
ここで、未知の〈Ｒr 2 （Ｅ）〉(3) を、上記の結果により既知となっった〈Ｒr 2 （Ｅ）〉(2) と摂動項Δr 2(3)を用いて、
〈Ｒr 2 （Ｅ）〉(3) ＝〈Ｒr 2 （Ｅ）〉(2) ＋Δr 2(3) ・・・（９９）
とおいて、上記の式（９８）を整理すると、下記の式（１００）となる。
【数７２】
【００７６】
したがって、摂動項Δr 2(3)は、下記の式（１０１）となる。
【数７３】
但し、式（１０１）におけるζr (3) （Ｅ）は下記の式（１０２）で表される。
【数７４】
【００７７】
以上のパラメータを用いて、２次の摂動モデルによるΔＲｐ(2) ，ΔＲpt(2) を、Ｌｉｎｄｈａｒｔが１次の摂動モデルで利用したように、
〈ΔＲp 2 （Ｅ）〉＝〈Ｒp 2 （Ｅ）〉＋〈Ｒp （Ｅ）〉2
＝｛〈Ｒc 2 （Ｅ）〉＋〈Ｒr 2 （Ｅ）〉｝／３−〈Ｒp （Ｅ）〉2 ・・・（１０３）
〈ＲT 2 （Ｅ）〉 ＝２｛〈Ｒc 2 （Ｅ）〉−〈Ｒr 2 （Ｅ）〉｝／３
・・・（１０４）
〈ΔＲpt2 （Ｅ）〉＝〈ＲT 2 （Ｅ）〉／２ ・・・（１０５）
で評価する。
但し、２次の場合には、１次と異なり、２次摂動モデルであることを意識して表記すると、上記の式（１０３）の第１式は、
〔〈ΔＲp 2 （Ｅ）〉(1) ＋Δp 2 (2) 〕
＝〔〈Ｒp 2 （Ｅ）〉(1) ＋Δp 2 (2) （Ｅ）〕
−〔〈Ｒp （Ｅ）〉(1) ＋Δp (2) （Ｅ）〕2 ・・・（１０６）
となる。
また、式（１０４）及び式（１０５）も同様である。
【００７８】
次に、３次の摂動モデルによるΔＲｐ(3) ，ΔＲpt(3) についても、Ｌｉｎｄｈａｒｔが１次の摂動モデルで利用したように、
〈ΔＲp 2 （Ｅ）〉＝〈Ｒp 2 （Ｅ）〉＋〈Ｒp （Ｅ）〉2
＝｛〈Ｒc 2 （Ｅ）〉＋〈Ｒr 2 （Ｅ）〉｝／３−〈Ｒp （Ｅ）〉2 ・・・（１０７）
〈ＲT 2 （Ｅ）〉 ＝２｛〈Ｒc 2 （Ｅ）〉−〈Ｒr 2 （Ｅ）〉｝／３
・・・（１０８）
〈ΔＲpt2 （Ｅ）〉＝〈ＲT 2 （Ｅ）〉／２ ・・・（１０９）
で評価する。
但し、３次の場合にも、１次と異なり、３次摂動モデルであることを意識して表記すると、上記の式（１０７）の第１式は、
〔〈ΔＲp 2 （Ｅ）〉(2) ＋Δp 2 (3) 〕
＝〔〈Ｒp 2 （Ｅ）〉(2) ＋Δp 2 (3) （Ｅ）〕
−〔〈Ｒp （Ｅ）〉(2) ＋Δp (3) （Ｅ）〕2 ・・・（１１０）
となる。
また、式（１０８）及び式（１０９）も同様である。
【００７９】
次に、３次のモーメントであるγを求めるために、式（３１）において、ｍ＝３とおくと、下記の式（１１１）が得られる。
【数７５】
【００８０】
また、上記のクロスタームに関する式（７３）において、ｍ＝１，ｋ＝２を代入すると、下記の式（１１２）が得られる。
【数７６】
【００８１】
ここで、下記の式（１１３）及び式（１１４）で表される変数を導入することによって、式（１１１）及び式（１１２）を整理することを試みる。
【数７７】
この変数の導入により、下記の式（１１５）及び式（１１６）のように、それぞれ変数が〈Ｒc 3 （Ｅ）〉及び〈Ｒr 3 （Ｅ）〉のみの独立な微分方程式を得ることができる。
【数７８】
【００８２】
この場合の変数〈Ｒc 3 （Ｅ）〉及び〈Ｒr 3 （Ｅ）〉は目的とするモーメントと、下記の式（１１７）及び式（１１８）によって関連付けられる。
【数７９】
【００８３】
一方、３次のモーメントμ3 は、Ｎ（ｘ）を規格化した濃度分布とすると、定義により、
μ3 ＝∫（ｘ−〈Ｒp （Ｅ）〉）3 Ｎ（ｘ）ｄｘ ・・・（１１９）
で表される。
したがって、３次のモーメントμ3 は、
μ3 ＝∫（ｘ−〈Ｒp （Ｅ）〉）3 Ｎ（ｘ）ｄｘ
＝∫（ｘ3 −３ｘ2 〈Ｒp （Ｅ）〉＋９ｘ〈Ｒp （Ｅ）〉2
−〈Ｒp （Ｅ）〉3 ）Ｎ（ｘ）ｄｘ
＝〈Ｒp 3 （Ｅ）〉−３〈Ｒp 2 （Ｅ）〉〈Ｒp （Ｅ）〉＋２〈Ｒp （Ｅ）〉3 ・・・（１２０）
となる。
【００８４】
したがって、
γ＝｛∫（ｘ−Ｒp （Ｅ））3 Ｎ（ｘ）ｄｘ｝／ΔＲp 3 ・・・（１２１）
で定義される分布の非対称性を表す３次のモーメントγは、下記の式（１２２）で評価される。
【数８０】
【００８５】
上記の微分方程式の解を下記の式（１２３）及び式（１２４）と置くと、
【数８１】
摂動項Δc 3(i)についての一般解が、式（１２５）として得られる。
【数８２】
但し、式（１２５）におけるζc 3(i)は、それぞれ、下記の式（１２６）乃至式（１２８）で表される。
【数８３】
【００８６】
また、他方の変数〈Ｒr 3 （Ｅ）〉に関しても、摂動項Δr 3(i)についての一般解が、式（１２９）として得られる。
【数８４】
但し、式（１２９）におけるζr 3(i)は、それぞれ、下記の式（１３０）乃至式（１３２）で表される。
【数８５】
【００８７】
次に、横方向の広がりΔＲptの深さ依存性について検討する。
横方向の広がりΔＲptは深さ依存性を有すると言われており、その依存性は、ΔＲpt0 をｘ＝Ｒp におけるΔＲptとすると、
ΔＲpt（ｘ）＝ΔＲpt0 ＋ｍ（ｘ−Ｒp ） ・・・（１３３）
で表現される。
【００８８】
そこで、〈ＲpT23（Ｅ）〉からｍの評価を試みる。
まず、ＲT 2 とΔＲptは、上述のように、
ＲT 2 ＝Δｘ2 ＋Δｙ2
＝２ΔＲpt2
＝２〔ΔＲpt0 ＋ｍ（ｘ−Ｒp ）〕2
≒２ΔＲpt02＋４ｍΔＲpt0 （ｘ−Ｒp ） ・・・（１３４）
で関連付けられる。
但し、最後の近似においては、ｘがＲp から遠い地点では濃度が急に減少するから、Ｒp 近傍のｍが重要であるため、ｍに関連する項は微小であると仮定している。
【００８９】
一方、〈ＲpT23（Ｅ）〉はｍと下記の式（１３５）の関係がある。
【数８６】
したがって、ｍは、
【数８７】
で表現することができる。
【００９０】
一方、Ｍｏｎｔｅ Ｃａｒｌｏ シミュレータによってイオン注入分布を発生させる場合には、毎回の衝突後に上記の式（３）の積分を実行するのではなく、衝突パラメータｂと散乱角度Φの関係式を数値的に解いて、その結果をグラフを用いてフィッティングにより解析的な近似式を導く。
【００９１】
この近似に際しては、一定高角度の項（π）と低角度の項で表しても良いし、一定高角度の項（π）、低角度の項及び中角度の項の３項で表しても良い。
この場合、低角度の項及び中角度の項を規格化衝突パラメータηを含む項と規格化エネルギーεの積で表す。
【００９２】
また、一定高角度の項（π）、低角度の項及び中角度の項をそれぞれ１より大きな数の累乗にして足し合わせたのち、その逆数の累乗をとって散乱角度Φを近似する。
【実施例１】
【００９３】
以上を前提として、ここで、図１乃至図４を参照して、本発明の実施例１のイオン注入分布発生方法を説明する。
図１は、本発明の実施例１のイオン注入分布発生方法のフローチャートであり、まず、ａ．基板種、注入不純物種、注入エネルギー、及び、ドーズ量からなる注入条件を入力す る。次いで、上記の式７５乃至式１０６で説明したように、
ｂ．２次の摂動モデルを用いて１次のモーメントである飛程の射影Ｒp と、２次のモーメ ントである射影Ｒp の偏差ΔＲp 及び横方向の偏差ΔＲptを求める。次いで、
ｃ．求めたＲp とΔＲp から１次元分布のガウス分布を発生させるとともに、Ｒp 、ΔＲ p 及びΔＲptから２次元のガウス分布を発生させる。
【００９４】
図２は、本発明の実施例１のイオン注入分布発生方法により求めたＲp とΔＲp を他のシミュレーションによって得たＲp とΔＲp と比較したものである。
ここでは、Ｍｏｎｔｅ Ｃａｒｌｏ法による結果とともに、ＬＳＳ理論による１次の摂動モデルによる結果、及び、上述のＧｉｂｂｏｎｓによる結果を１．５次の摂動モデルとして示している。
【００９５】
図から明らかなように、１次のモーメントである飛程の射影Ｒp については、全てのシミュレーション結果が良好な一致を見せている。
しかし、２次のモーメントである射影Ｒp の偏差ΔＲp については、ＬＳＳ理論による１次の摂動モデルによる結果に比べて、本発明の実施例１の２次の摂動モデルによる結果はＭｏｎｔｅ Ｃａｒｌｏ法による結果と非常に良好な一致が見られる。
なお、１．５次の摂動モデルは、射影Ｒp の偏差ΔＲp については、積分方程式の分母に発散項を含んでいるため、計算不能となっている。
【００９６】
図３は、本発明の実施例１のイオン注入分布発生方法により求めた１次元分布であり、ここでは、Ｓｉ基板にＰ（リン）イオンを１０ｋｅＶで１×１０15ｃｍ-2注入した場合の結果を示している。
図から明らかなように、綺麗なガウス分布を示している。
【００９７】
図４は、本発明の実施例１のイオン注入分布発生方法により求めた２次元分布であり、ここでも、Ｓｉ基板にＰ（リン）イオンを１０ｋｅＶで１×１０15ｃｍ-2注入した場合の結果を示している。
なお、この場合には、イオン注入マスクはイオンを全く透過しないものとしてシミュレーションしている。
【００９８】
このように、本発明の実施例１においては、２次の摂動モデルを用いることにより、簡単な計算により、Ｍｏｎｔｅ Ｃａｒｌｏ法に比べて極めて短時間で射影Ｒp 、偏差ΔＲp 、横方向の偏差ΔＲptをＭｏｎｔｅ Ｃａｒｌｏ法に近い精度で求めることができる。
【００９９】
なお、この実施例１においては、２次の摂動モデルであるので、２次のモーメントまでしか精度良く求めることができないので、発生できるイオン注入分布はガウス分布となる。
しかし、実際にＬＳＳプロセスシミュレーションで標準的に利用されているＰｅａｒｓｏｎ分布は、γやβ等のより高次のモーメント情報が必要である。
【実施例２】
【０１００】
ここで、３次及び４次のモーメントを求めて、より実測結果に近いＰｅａｒｓｏｎ分布を発生させるために、３次の摂動モデルを用いた本発明の実施例２のイオン注入分布発生方法を説明する。
図５は、本発明の実施例２のイオン注入分布発生方法のフローチャートであり、まず、ａ．基板種、注入不純物種、注入エネルギー、及び、ドーズ量からなる注入条件を入力す る。次いで、上記の式７９乃至式１３５で説明したように、
ｂ．３次の摂動モデルを用いて１次のモーメントである飛程の射影Ｒp と、２次のモーメ ントである射影Ｒp の偏差ΔＲp 及び横方向の偏差ΔＲptを深さ依存性を有するよう にを求める。また、３次のモーメントであるγを求めるとともに、４次のモーメント であるβを経験的に求める。次いで、
ｃ．求めたＲp 、ΔＲp 、γ及びβからから１次元分布を発生させるとともに、Ｒp 、Δ Ｒp 、γ、β及び深さ依存性ΔＲptから２次元分布を発生させる。
【０１０１】
ここでは、４次のモーメントであるβの求め方を説明する。
βも他のモーメントと同様にＬＳＳ積分方程式を解くことによって原理的には求めることができるが、３次のモーメントγの場合よりも更に複雑になり、なおかつ数値的に不安定である。
そこで、イオン注入分布はそれほど特異な変形をしないため、βとγの間にはユニバーサルが関係があることをＭｏｎｔｅ Ｃａｒｌｏで示し、それを経験式で表現したものを利用する。
【０１０２】
図６は、Ｍｏｎｔｅ Ｃａｒｌｏによるシミュレーション結果をβとγとの関係として示したものである。
なお、ここでは、Ｓｉ基板に対してＢ、Ａｓ、及び、Ｐをそれぞれ、１０ｋｅＶ、２０ｋｅＶ、４０ｋｅＶ、８０ｋｅＶ、１６０ｋｅＶ、及び、３２０ｋｅＶでイオン注入した場合の結果を示している。
【０１０３】
この相関関係は、
β＝２．６６１＋１．８５２γ2 ・・・（１３７）
と表現される。
この近似式を用いて、３次の摂動モデルにより求めたγからβを求めて、上述の１次元分布及び２次元分布を発生させれば良い。
【０１０４】
但し、この関係式（１３７）は、図に示すようによく利用されるＰｅａｒｓｏｎ ＩＶの領域になく、Ｐｅａｒｓｏｎ Ｉ、ＩＩ、ＩＩＩ、ＶＩの領域を横切っている。
そこで、図にβD2として示すように、βを５％増加させて強制的にＰｅａｒｓｏｎ ＩＶの下記の式（１３７）で示される臨界値βD2に移動させて用いる。
βD2＝〔３（１３γ2 ＋１６）＋６（γ2 ＋４）3/2 〕／（３２−γ2 ）
・・・（１３８）
【実施例３】
【０１０５】
次に、衝突パラメータｂと散乱角度Φの関係式を解析的に表現して計算コストを低減する本発明の実施例３のイオン注入分布発生方法を説明する
上記の式（１）において、伝達エネルギーＴ2fが、ある程度小さくなれば、例えば、１０ｅＶ程度になれば、それ以上大きなｂは考える必要がない。
つまり核との相互作用はないと考えることができる。
【０１０６】
このような伝達エネルギーをＴ2fc とすると、この場合の散乱角度Φc も小さいと考えることができるので、ｓｉｎ2 （Φ／２）はΦc 2 ／４で近似することができる。
したがって、上記の式（１）は、
Ｔ2fc ／Ｔ1i ＝〔４Ｍ2 Ｍ1 ／（Ｍ1 ＋Ｍ2 ）2 〕×Φc 2 ／４ ．．（１３９）８なる。これより散乱角度Φc は
Φc ＝｛〔（Ｍ1 ＋Ｍ2 ）2 ／Ｍ2 Ｍ1 〕×（Ｔ2fc ／Ｔ1i）｝1/2 ・・（１４０）と求まる。
ここで、散乱角度Φc に対応する衝突パラメータをｂmax とすると、散乱角度Φc は下記の式（１４１）で表される。
【数８８】
【０１０７】
実際は衝突のたびにこの積分からｂmax を求めるのは計算コストがかかるため、予め式（３）からｂとΦとの関係式を数値的に解いて、その結果をフィッティングにより解析的に表現する式を用いる。
なお、フィッティングに際しては、上記の式（３）と同様なユニバーサルな変数で示した下記の式（１４２）を用いる。
【数８９】
但し、ここでは、 η＝ｂmax ／ａU ・・・（１４３）
であり、ｂmax に対応するρがρmin である。
【０１０８】
図７は、散乱角度Φのエネルギーε依存性の説明図であり、上記の式（１４２）を数値的に解いた結果を実線で示している。
ここで、各ηについて、フィッティングの結果、
Φc ＝〔（１／π）＋１／（ａ（η）ε-0.98 ）〕-1 ・・・（１４４）
で合わせ込んだ結果を破線で示している。
【０１０９】
図に示すように、角度が大きくなって飽和に向かう領域でフィッティングはうまくいっていない。しかし、実際に利用するのは数度の角度、つまり１０-2ｒａｄｉａｎ の領域であるから問題ない。また、この領域は簡便に、
Φc ＝ａ（η）ε-0.98 ・・・（１４５）
と表現される。ただし、ａ（η）は下記の式（１４６）で表される。
【数９０】
【０１１０】
しかし、実際に利用するにはηをａであらわしたη（ａ）が便利であるので、ａ−ηの関係式からフィッティングする。
図８は、上記の式（１４６）におけるａのη依存性を説明図であり、これを逆のηのａ依存性として示したのが図９である。
この図９からフィッティングにより求めたη（ａ）が下記の式（１４７）で表される。
【数９１】
これから、
ａ＝Φc ／ε-0.98 ・・・（１４８）
の場合のηが、上記の式（１４６）から求まる。
さらに、ηが求まると、上記の式（１４２）からｂmax が求まる。
【０１１１】
核との相互作用の衝突断面積σn は
σn ＝πｂmax 2 ・・・（１４９）
と考えることができる。
ここで、基板原子の濃度をＮとすると：
ＬＮσn ＝１ ・・・（１５０）
となる距離Ｌだけ進むと一回衝突すると考えることができる。
【０１１２】
そこで、エネルギーＥのイオンが距離Ｌの間は核との相互作用なしに進み、その後に核と相互作用する、と現象を記述する。Ｌ進む間は電子とのみ相互作用すると考える。
即ち、Ｌ進む間は電子との相互作用で
ΔＥc ＝Ｓe （Ｅ）Ｌ ・・・（１５１）
のエネルギーを失う。
Ｌ進んだ後の核との相互作用は衝突パラメータｂを０からｂmax までの円上の点をランダムに選ぶ。つまり
ｂ＝ｂmax 〔Ｒａｎｄ（ｎ）〕1/2 ・・・（１５２）
とする。
但し、Ｒａｎｄ（ｎ）はｎが０から１の間の乱数である。
【０１１３】
Ｌは原子間の平均距離Ｌmin （＝１／Ｎ1/3 ）より小さくなった場合は、
Ｌ＝Ｌmin ・・・（１５３）
としなければならない。この場合、衝突と関連する面積は矩形で、 σn ＝Ｌmin 2 ・・・（１５４）
となる。
これをこれまで扱ってきた衝突パラメータに焼きなおすと、この場合のｂmax は、上記の式（１４９）から、
ｂmax ＝Ｌmin ／π1/2 ・・・（１５５）
となる。
したがって、このｂmax を使って上記の式（１５２）にしたがって乱数を発生させてｂを定めて行けば良く、注入エネルギーが上記の式（１４１）が精度良く成立する予め定めた最低臨界エネルギーに達するまで計算を行う。
【０１１４】
図１０は、１０個の粒子の軌跡の各Ｌを示したもので、ここでは、Ｌmin で規格化しており、ＬはＬmin の１０倍程度から徐々に小さくなってくる。
したがって、常にＬmin を用いた場合より計算回数が著しく少なくなり、計算時間が短縮される。
なお、図においては１０個の粒子の各Ｌを各々Ｌ１〜Ｌ１０で示している。
【０１１５】
図１１は、イオン注入分布の説明図であり、ここでは、Ｂイオンを４０ｋｅＶの加速エネルギーで、１×１０15ｃｍ-2注入した場合のイオン注入分布を示している。
図から明らかなように、本発明の実施例３の可変Ｌで求めたイオン注入分布は、常に、Ｌmin を用いた場合と計算結果はほとんど変わらず、且つ、ＳＩＭＳによる実測結果を良好に表している。
【０１１６】
ここで、複数の原子から基板が構成されている場合の取り扱いを考える。
原子ｉがｎi 個あるとする。それぞれの原子の臨界角度Φciは
Φci＝｛〔（Ｍ1 ＋Ｍ2i）2 ／Ｍ1 Ｍ2i〕×（Ｔ2fc ／Ｔ1i）｝1/2 ・・（１５６）と求まり、これから衝突断面積σniが計算でき、トータルの衝突断面積σn は、下記の式（１５６）で表される。
【数９２】
【０１１７】
ここで、基板原子の濃度をＮとすると、
ＬＮσn ＝１ ・・・（１５８）
となる距離Ｌだけ進むと一回衝突すると考えることができる。
【０１１８】
その衝突で原子ｉと衝突する確率は、下記の式（１５９）で表される。
【数９３】
そこで、その確率に対応させて乱数を発生させて原子ｉを選択する。選択した後は構成原子が一種類の場合と同じである。
また、Ｌが、下記の式（１６０）より小さくなった場合には、Ｌ＝Ｌmin としなければならない。
【数９４】
【０１１９】
また、この場合は衝突断面積は原子密度で決定されるから、その衝突で原子ｉと衝突する確率は２種類の場合と同様と考えると、ｂmaxiは、
ｂmaxi＝（ｒσni／π）1/2 ・・・（１６１）
となる。但し、ｒは下記の式（１６２）で表される。
【数９５】
【０１２０】
このように、本発明の実施例３においては、散乱角度が小さい場合に場合に、積分からではなく、フィッティングによって求めた近似式を用い、且つ、可変Ｌを用いて解析的にイオン注入分布を求めているので、計算精度を落とすことなく計算回数を著しく少なくすることができる。
【実施例４】
【０１２１】
次に、散乱角度が大きい場合にも有効であるように、遷移領域での精度を向上させるために、式（１４４）を拡張した本発明の実施例４のイオン注入分布発生方法を説明する。
図１２は、散乱角度Φのエネルギーε依存性の説明図であり、上記の式（３）を数値的に解いた結果を実線で示している。
また、図１３は、図１２における１〜３ラジアンの範囲を拡大した図である。
ここで、各ηについて、フィッティングの結果、
Φ＝１／｛（１／π）3/2 ＋〔１／（ａ（η）ε-0.98 ）〕3/2
＋〔１／ｇ（η）ε-0.2〕3/2 〕2/3 ・・・（１６３）
で合わせ込んだ結果を破線で示している。
【０１２２】
なお、ここで、中間の散乱角度に対応する項ｇ（η）は、
ｇ（η）＝〔１／（４．０η-0.46 ）＋１／（９．０η-1.3）〕-1
＝４．０η-0.46 ／（１＋０．４４η0.84） ・・・（１６４）
で表され、図１４に示すη依存性が見られる。
【０１２３】
図１２及び図１３から明らかなように、低角度から高角度まで精度良く積分による数値計算の結果を再現している。
なお、上記の式（１６３）における（１／π）3/2 が高角度の項に相当する。
【図面の簡単な説明】
【０１２４】
【図１】本発明の実施例１のイオン注入分布発生方法のフローチャートである。
【図２】本発明の実施例１のイオン注入分布発生方法により求めたＲp とΔＲp の説明図である。
【図３】本発明の実施例１のイオン注入分布発生方法により求めた１次元分布である。
【図４】本発明の実施例１のイオン注入分布発生方法により求めた２次元分布である。
【図５】本発明の実施例２のイオン注入分布発生方法のフローチャートである。
【図６】Ｍｏｎｔｅ Ｃａｒｌｏによるシミュレーション結果をβとγとの関係として示した図である。
【図７】散乱角度Φのエネルギーε依存性の説明図である。
【図８】式（１４６）におけるａのη依存性を説明図である。
【図９】ηのａ依存性を説明図である。
【図１０】各軌跡のＬの説明図である。
【図１１】イオン注入分布の説明図である。
【図１２】散乱角度Φのエネルギーε依存性の説明図である。
【図１３】図１２における１〜３ラジアンの範囲の拡大図である。
【図１４】ｇのη依存性の説明図である。
【図１５】計算モデルである。
【図１６】イオンの飛程Ｒの模式図である。
【図１７】Ｒp （Ｅ，ｃｏｓφ），ＲT （Ｅ，ｃｏｓφ），Ｒp （Ｅ），ＲT （Ｅ）の幾何学的関係の説明図である。
【技術分野】
【０００１】
本発明はイオン注入分布発生方法及びシミュレータに関し、Ｍｏｎｔｅ Ｃａｒｌｏ ｓｉｍｕｌａｔｏｒとほぼ同じ精度で、より短時間でイオン注入分布を発生させるための手法に関するものである。
【背景技術】
【０００２】
シリコン集積回路装置において、シリコン基板への不純物の導入はイオン注入で行われるのが一般的である。
このようなシリコン集積回路装置のプロセス構築に際しては、必要な素子構造を得るためのイオン注入条件を決定する必要があるが、このようなイオン注入条件をシミュレーションにより決定することが行われている。
非晶質層へのイオン注入分布を理論的に予想する手段としてＭｏｎｔｅ Ｃａｒｌｏ がある。これは、入射イオンと基板との相互作用を、核阻止能及び電子阻止能の物理に基づいて、入射イオンの軌跡を追跡していくものである。
【０００３】
この理論は、任意のイオンを任意の基板にイオン注入した場合の一般的な場合にも有効であり、電子阻止能をチューニングすればその精度をさらに向上させることができる。
図１５は、計算モデルであり、質量数Ｍ1 ，原子番号Ｚ1 ，エネルギーＴ1i（速度ｖ1i）のイオンが、基板を構成する質量数Ｍ2 ，原子番号Ｚ2 の原子と相互作用して伝達するエネルギーＴ2fは、散乱角度をΦ、相互作用後のイオンの速度をν1i、基板原子の速度をν2iとすると、
Ｔ2f／Ｔ1i＝２Ｍ2 ν2i2 ｓｉｎ2 （Φ／２）／〔（１／２）Ｍ1 ν1i2 〕
＝（４Ｍ2 ／Ｍ1 ）｛〔Ｍ1 ｖ1i／（Ｍ1 ＋Ｍ2 ）〕2 ／ν1i2 ｝
×ｓｉｎ2 （Φ／２）
＝〔４Ｍ2 Ｍ1 ／（Ｍ1 ＋Ｍ2 ）2 〕ｓｉｎ2 （Φ／２）・・・（１）
と表現される。
【０００４】
ここで、イオンと基板原子の距離をｒ、衝突パラメータをｂ、ポテンシャルエネルギーをＶ（ｒ）とすると、伝達エネルギーＴ2fは、下記の式（２）として求まる。
【数１】
つまり、注入されたイオンは、核との相互作用により、
ΔＥn ＝Ｔ2f
のエネルギーを失う。
【０００５】
また、相互作用に伴う散乱角度Φは、下記の式（３）で表される。
【数２】
ここで、半径ｂの円周上の位置、即ち角度θは、Ｒａｎｄ（ｎ）をｎが０から１の間の乱数とすると、
θ＝２πＲａｎｄ（ｎ）
の関係から求まる。
【０００６】
これにより、衝突後のエネルギーと方向が決まる。
次に、新たに注入エネルギーで同様の計算を繰り返し、イオンの軌跡をトレースしていけば良い。
しかし、実際の衝突のたびに、上記の式（３）で表される積分を毎回実行するのは計算コストが膨大になるため、ＺｉｅｇｌｅｒはＭａｇｉｃ ｆｏｒｍｕｌａを提案し、この積分を簡単なプロトコルで解いている。
また、シミュレータＴＳＵＰＲＭでは、計算結果をテーブル化しその内挿でこの値を評価している。
【０００７】
なお、上記の式（３）における各パラメータを規格化したユニバーサルな変数で表すと、下記の式（４）で表される。
【数３】
但し、ａB をボーア径とすると、
ρ＝ｒ／ａU ，
η＝ｂ／ａU ，
ａU ＝０．８８５４ａB ／（Ｚ10.23 ＋Ｚ20.23 ）
である。
また、エネルギーは、
ε＝ＥC ／（Ｚ1 Ｚ2 ｅ2 ／ａU ）
Ｅc ＝（１／２）Ｍc ｖ1i2
１／Ｍc ＝１／Ｍ1 ＋１／Ｍ2
である。
また、ポテンシャルエネルギーＶ（ｒ）として、下記のＺｉｅｇｌｅｒ−Ｌｉｔｍａｒｋ−Ｂｉｅｒｓａｋ（ＺＬＢ）のポテンシャルエネルギーを用いる（例えば、非特許文献１参照）。
Ｖ（ｒ）＝（ｅ2 Ｚ1 Ｚ2 ／ｒ）ｆ（ρ）
但し、
【数４】
【０００８】
しかしながら、上述の計算コスト低減の手法を導入しても、Ｍｏｎｔｅ Ｃａｒｌｏシミュレーションは、粒子の各軌跡を追うため、統計誤差を減らすためには数万個以上の計算をする必要があり、時間がかかるという問題がある。
【０００９】
そこで、電子阻止能Ｓe 、核阻止能Ｓn が与えられた時、イオン注入分布が従うべき積分方程式が、Ｌｉｎｄｈａｒｔ，Ｓｃｈａｒｆ，Ｓｃｈｉｏｔｔによって提案され（例えば、非特許文献２参照）、このモデルはＬＳＳ理論と呼ばれている。
【００１０】
この理論では、粒子の軌跡を追跡することなしに注入条件が決まれば、積分方程式を解くことによって分布の任意の次数の分布モーメントのエネルギー依存性までを即座に計算できる。
【００１１】
この場合、積分方程式を展開し、微分方程式に還元し、それを解いていく必要があり、そのステップで近似が入る。
即ち、伝達エネルギーが初期エネルギーに比べて小さいと仮定してモーメントをエネルギーに関してＴａｙｌｏｒ展開するが、この展開する次数が解くべき微分方程式の階数となる。
【００１２】
また、散乱角度もエネルギーに関してＴａｙｌｏｒ展開するが、これは解くべき方程式の係数と関連し階数とは無関係である。解析モデルは一階の線形微分方程式の場合のみ任意の係数で一般的に得られる。
このＬＳＳ理論ではモーメントに関しては１次、散乱角度に関しては２次までＴａｙｌｏｒ展開し、線形１階の微分方程式を解きＲp 及びΔＲp の解析モデル式を導出している。
【００１３】
ここで、図１６を参照して、ＬＳＳ理論における注入イオンの飛程Ｒに関する積分方程式を説明する。
図１６はイオンの飛程Ｒの模式図であり、エネルギーＥで注入されたイオンは基板原子と相互作用しながらエネルギーを失い方向を変えながら図に示すように進んでいき、全エネルギーを失い基板中で静止する。
【００１４】
この時、イオンが進んだ軌跡の線積分を飛程Ｒとし、飛程Ｒの注入方向への射影をＲp 、横方向への広がりをＲT とする。
また、横方向への広がりのｘ成分をΔｘ、ｙ成分をΔｙとする。
なお、横方向への広がりは、図面及び以降の積分方程式或いは微分方程式においてはＲに「垂直記号」のサフィックスを付けた記号で表すが、明細書本文中では、明細書作成の都合上、「ＲT 」で表す。
【００１５】
表面に垂直にエネルギーＥでイオン注入されたイオンがＲとＲ＋ｄＲの間で止まる確率をＰ（Ｅ，Ｒ）とする。この場合、
【数５】
である。
また、Ｒに関してｍ次のモーメントを
【数６】
と定義する。
【００１６】
イオンが表面からｄＲ進む間に原子および電子と相互作用する確率は、
ＮｄＲ∫ｄσn ＋ＮｄＲ∫ｄσe ・・・（８）
であり、そのそれぞれの相互作用でイオンのエネルギーは、
Ｅ−Ｔn ，Ｅ−Ｔe
に減少すると仮定する。
ここで、ｄσn ，ｄσe は核阻止能及び電子阻止能と関連する微分断面積である。
【００１７】
よって、このイオンが原子および電子と相互作用し、エネルギーを失って飛程Ｒに止まる確率は、
【数７】
で与えられる。
また、ｄＲ進む間に衝突しない確率は、
１−（ＮｄＲ∫ｄσn ＋ＮｄＲ∫ｄσe ） ・・・（１０）
である。この間にエネルギーは失われないから、イオンがＲに止まる確率は、
［１−（ＮｄＲ∫ｄσn ＋ＮｄＲ∫ｄσe ）］Ｐ（Ｅ，Ｒ−ｄＲ） ・・・（１１）
となる。よって、
【数８】
となる。
【００１８】
ここで、ｄＲを無限小に持っていく極限を考えると、
【数９】
となる。両辺にＲm をかけてＲについて０から∞まで積分すると、
【数１０】
となる。
【００１９】
式（１４）の左辺を部分積分し、また右辺の積分の順序を入れ換えると、
【数１１】
式（１５）における左辺の第１項は０である。また
【数１２】
と表現するが、これはＲm の期待値に相当する。
【００２０】
ここで、式（１６）を用いると、式（１５）は、
【数１３】
となり、これが飛程Ｒに関するｍ次の積分方程式である。
ここで、飛程Ｒに関しては１次に関してのみ解析する。式（１７）においてｍ＝１とおいて
【数１４】
を得る。これが、飛程Ｒの従うべき積分方程式である。
【００２１】
次に、式１７をＴn ，Ｔe ≪Ｅの近似のもとで１次の項までを残し解いていく。この近似は伝達エネルギーが小さいということを仮定をしている。
この場合、核との相互作用に関しては広角散乱を無視する近似に相当し、電子との相互作用に関しては、エネルギーが高いほどＴe は大きくなる。このため、比較的エネルギーの低い場合に近似の精度は良くなってくる。
この解を１次のものであることを意識して〈Ｒ（Ｅ）〉(1) と表現する。
【００２２】
次に、式（１８）の両辺をＮで割るとともにＴａｙｌｏｒ展開することにより、
【数１５】
となる。ここで、
Ｓn ＝∫Ｔn ｄσn ，Ｓe ＝∫Ｔe ｄσe ・・・（２０）
を利用している。式（１９）より、良く知られている、
【数１６】
が導出される。
【００２３】
次に、図１７を参照して、Ｒp （Ｅ，ｃｏｓφ），ＲT （Ｅ，ｃｏｓφ），Ｒp （Ｅ），ＲT （Ｅ）の幾何学的関係を説明する。
図１７は、Ｒp （Ｅ，ｃｏｓφ），ＲT （Ｅ，ｃｏｓφ），Ｒp （Ｅ），ＲT （Ｅ）の幾何学的関係の説明図であり、エネルギーＥ、角度φで入射したイオンが基板内のＡ点に静止した状況を模式図的に示している。
入射角度０の軸に垂直な面上での横方向の広がりをＲT （Ｅ，ｃｏｓφ）とすると、入射方向に垂直な軸に対する射影Ｒp （Ｅ）が、入射方向に垂直な面に対する横方向広がりがＲT （Ｅ）と考えることができる。よって、
【数１７】
【００２４】
以上の幾何学的関係を基にして、次に、Ｒp ，ＲT ，Ｒp ＲT の従うべき積分方程式を導出していく。
エネルギーＥで、角φで散乱されたイオンの射影飛程がＲp とＲp ＋ｄＲp の間で止まる確率をＰp （Ｅ，Ｒp ，ｃｏｓφ）とする。
φ＝０の場合の確率をＰp （Ｅ，Ｒp ）とする。また、Ｒp に関してｍ次のモーメントを、
【数１８】
と定義する。
【００２５】
入射イオンがｄＲp 進む間に原子および電子と相互作用する確率は、飛程Ｒの場合と同様に、
ＮｄＲp ∫ｄσn ＋ＮｄＲp ∫ｄσe ・・・（２５）
であり、この相互作用でのエネルギーは、
Ｅ−Ｔn ，Ｅ−Ｔe
に減少する。
【００２６】
この場合、電子との相互作用による散乱角は０とみなすので、このイオンが原子および電子と相互作用し、エネルギーＴn ，Ｔe を失ってＲp に止まる確率は、
【数１９】
で与えられる。
また、ｄＲ進む間に衝突しない確率は、
１−（ＮｄＲp ∫ｄσn ＋ＮｄＲp ∫ｄσe ） ・・・（２７）
である。この間にエネルギーは失われないから、イオンがＲp に止まる確率は、
［１−（ＮｄＲp ∫ｄσn ＋ＮｄＲp ∫ｄσe ）］Ｐ（Ｅ，Ｒp −ｄＲp ）
・・・（２８）
となる。よって、
【数２０】
となる。
【００２７】
ここで、ｄＲp を無限小にもっていく極限を考えると、
【数２１】
となり、両辺にＲm p をかけてＲp について０から∞まで積分すると、
【数２２】
とｍ次の〈Ｒm p （Ｅ）〉に関する積分方程式が導出される。
【００２８】
次に、分布の横方向への広がりを表すＲT について考える。
エネルギーＥ、入射角φで注入されたイオンがＲT とＲT ＋ｄＲT の間に止まる確率をＰT （Ｅ，ＲT ，ｃｏｓφ）とする。
この時、入射角φのイオンがｄＲT 進む間に原子および電子と相互作用する確率は、
【数２３】
である。よって、このイオンが相互作用しエネルギーＴn ，Ｔe を失ってＲT に止まる確率は、
【数２４】
で与えられる。
【００２９】
この場合も電子との相互作用では角度の変化はなく、さらにこの場合はＲT も変化しないことを仮定している。また、ｄＲT 進む間に衝突しない確率は、
【数２５】
である。この間にエネルギーは失われず、角度も変化しないためＲT も変化しない。よって、イオンがＲT に止まる確率は、
【数２６】
となる。
【００３０】
以上より、
【数２７】
となる。
【００３１】
ここで、ｄＲT を無限小にもっていく極限を考えると、
【数２８】
となり、両辺にＲm T をかけてＲT について０から∞まで積分すると、
【数２９】
と〈Ｒm T （Ｅ）〉に関する積分方程式が導出される。
【００３２】
まず、飛程の射影Ｒp を求めるために、式（３１）においてｍ＝１とすると下記の式（３９）となる。
【数３０】
ここで、１次のテーラー展開した解を〈Ｒp （Ｅ）〉(1) とおくと、下記の式（４０）となる。
【数３１】
【００３３】
これから、下記の式（４１）が得られる。
【数３２】
但し、α（Ｅ）及びα（Ｅ′）は、下記の式（４２）で表される。
【数３３】
【００３４】
次に、２次のモーメントΔＲp 及びΔＲptを求めるために、式（３１）及び式（３８）においてｍ＝２とすると下記の式（４３）及び式（４４）となる。
【数３４】
【００３５】
ここで、Ｌｉｎｄｈａｒｔは下記の式（４５）及び式（４６）で定義される変数〈Ｒ2 c （Ｅ）〉及び〈Ｒ2 r （Ｅ）〉を導入した。
【数３５】
【００３６】
〈ΔＲ2 p （Ｅ）〉と〈Ｒ2 T （Ｅ）〉は、変数〈Ｒ2 c （Ｅ）〉及び〈Ｒ2 r （Ｅ）〉と下記の式（４７）及び式（４８）で関連づけられる。
【数３６】
【００３７】
また、横方向の偏差〈ΔＲ2 pt（Ｅ）〉は〈ΔＸ2 〉と関連し、〈Ｒ2 T （Ｅ）〉そのものではない。したがって、
〈Ｒ2 T （Ｅ）〉＝〈ΔＸ2 〉＋〈ΔＸ2 〉＝２〈ΔＲ2 pt（Ｅ）〉・・・（４９）
となり、これより、
〈ΔＲ2 pt（Ｅ）〉＝〈Ｒ2 T （Ｅ）〉／２ ・・・（５０）
となる。
【００３８】
Ｌｉｎｄｈａｒｔはこれらの変数を導入することにより、下記の式（５１）及び式（５２）の独立な微分方程式を得ることに成功した。
【数３７】
【００３９】
また、Ｇｉｂｂｏｎｓは、このモデル式にエネルギーに依存しない核阻止能を代入し、Ｓｉ中のＢ，Ｎ，Ａｌ，Ｐ，Ｇａ，Ａｓ，Ｉｎ，ＳｂのＲp ，ΔＲp を評価し、いくつかの実験データと比較し、良い一致を得ている。
さらに、３次までのモーメントを数値的に評価することが検討されている。
【非特許文献１】Ｊ．Ｆ．Ｚｉｅｇｌｅｒ，Ｊ．Ｐ．Ｂｉｅｒｓａｃｋ，ａｎｄ Ｕ．Ｌｉｔｔｍａｒｋ，Ｔｈｅ ｓｔｏｐｐｉｎｇ ａｎｄ ｒａｎｇｅ ｏｆ ｉｏｎｓ ｉｎ ｓｏｌｉｄ，Ｐｅｒｇａｍｏｎ，１８８５
【非特許文献２】Ｊ．Ｌｉｎｄｈａｒｔ，Ｍ．Ｓｃｈａｒｆｆ，Ｈ．Ｓｃｈｉｏｔｔ，Ｍａｔ．Ｆｔｓ．Ｍｅｄｄ．Ｖｉｄ．Ｓｃｌｓｋ，ｖｏｌ．３３，ｐｐ．１−３９，１９６３
【発明の開示】
【発明が解決しようとする課題】
【００４０】
しかし、ＬＳＳ理論では、積分方程式を解く際に近似が必要で、これまでに提案されているモデルではＲp は問題ないが、ΔＲp は格段に精度が落ちるか、計算できないという問題がある。
【００４１】
また、上述のように、３次までのモーメントを数値的に評価することが検討されているが、これらの数値的に解析したデータは系統的に比較されていない。また、比較するうえでも、３次までのモーメントだけでは現在利用されているＰｅａｒｓｏｎ分布を利用することが出来ないという問題点もある。
【００４２】
即ち、核阻止能モデル、電子阻止能モデルの種類、積分方程式の展開の仕方によってＬＳＳ理論から導かれる結果はさまざまになる。ＬＳＳ理論で扱われている近似の精度は同じ核阻止能モデル、電子阻止能モデルを使ったＭｏｎｔｅ Ｃａｒｌｏ計算と比較することによって正確に評価される。
【００４３】
しかしながら、精度の高い核阻止能、電子阻止能を使ってＭｏｎｔｅ ＣａｒｌｏとＬＳＳ理論の比較がこれまでなされていない。さらに実際の分布は４次のモーメントまで利用して初めて精度よく表現される。
したがって、これまでのＬＳＳ理論では、現在利用されているＰｅａｒｓｏｎ関数と対応づけることはできない。
【００４４】
また、衝突パラメータｂと散乱角度Φの相関からイオン注入分布を求める場合も、Ｍｏｎｔｅ Ｃａｒｌｏは粒子の各軌跡を追うため、統計誤差を減らすためには数万個以上の計算をする必要があるため時間がかかる。
この計算を進めていく上で一番時間のかかるのが散乱角度を評価するための積分である。
【００４５】
また、衝突と次の衝突の間の距離Ｌをいくらにとるのかも計算時間に影響をおよぼす。
これまでは平均原子間距離Ｌmin が利用されており、この場合は相互作用が無視できるような場合でも散乱を計算してしまうため計算時間が膨大になるという問題がある。
【００４６】
したがって、本発明は、実験データを再現するＭｏｎｔｅ Ｃａｒｌｏ ｓｉｍｕｌａｔｏｒと同じ核阻止能、電子阻止能を用いたＬＳＳ理論の積分方程式から、２次、３次の摂動計算によってＲp ，ΔＲp ，ΔＲpt，γの解析モデルを導出して計算時間を大幅に短縮することを目的とする。
【００４７】
また、Ｍｏｎｔｅ Ｃａｒｌｏを用いる場合にも、散乱角度を評価するための積分を解析的に表現する近似式で置き換えて計算時間を大幅に短縮することを目的とする。
【課題を解決するための手段】
【００４８】
本発明の一観点によると、半導体に注入するイオンの飛程の射影Ｒp の２次の項まで考慮した射影〈Ｒp （Ｅ）〉(2) を、１次の項まで考慮した既知の射影を〈Ｒp （Ｅ）〉(1) とした時に、摂動項Δp (2) （Ｅ）を用いて、
〈Ｒp （Ｅ）〉(2) ＝〈Ｒp （Ｅ）〉(1) ＋Δp (2) （Ｅ）
とした近似式を用いた２次の摂動モデルを用いて求めるイオン注入発生方法が提供される。
【００４９】
また、別の観点によると、基板に注入したイオンの散乱角度Φをモンテ・カルロ法で求める際に、散乱角度Φを、飽和する一定高角度の項、低角度の項、中間角度の項の３つの項の和の逆数で解析式として表すイオン注入分布発生方法が提供される。
【００５０】
さらに、別の観点によると、上述の各種のイオン注入分布発生方法に関する計算式が組み込まれたシミュレータが提供される。
【発明の効果】
【００５１】
本発明によれば、精度の高い解析モデルを導出し、さらに高次のモーメントを高精度で予想できるようになった。これにより、現在標準的に用いられているＰｅａｒｓｏｎ分布のパラメータを初めて解析的に求めることが可能になる。
したがって、任意不純物を任意非晶質基板にイオン注入した分布のエネルギー依存性が瞬時に予想できるようになる。
また、この分布発生ソフトにより横方向の分布の深さ依存性も予想できる。
【００５２】
さらに、Ｍｏｎｔｅ Ｃａｒｌｏを用いる場合にも、散乱角度を評価するための積分を解析的に表現する近似式で置き換えてることによって、計算時間を大幅に短縮することが可能になる。
【発明を実施するための最良の形態】
【００５３】
まず、１次乃至３次の各阻止能について、計算しておく。
１次の核阻止能Ｓn （Ｅ）は、下記の式（５３）で与えられる。
【数３８】
また、ユニバーサルな核阻止能Ｓn （ε）は、下記の式（５４）で与えられる。
【数３９】
【００５４】
ここで、上述のＺＬＢポテンシャルを利用すると、ユニバーサルな核阻止能Ｓn （ε）は、下記の式（５５）で近似される。
【数４０】
【００５５】
この阻止能は高次のものにも拡張され、２次のモーメントに対応する拡張された核阻止能Ωn 2 は、下記の式（５６）で表される。
【数４１】
但し、ユニバーサルな拡張された核阻止能Ωn 2 （ε）は下記の式（５７）で表される。
【数４２】
これも、下記の式（５８）で近似的に表現される。
【数４３】
【００５６】
さらに、３次のモーメントに対応する拡張された核阻止能Λn 3 は、下記の式（５９）で表される。
【数４４】
但し、ユニバーサルな拡張された核阻止能Λn 3 （ε）は下記の式（６０）で表される。
【数４５】
これも、下記の式（６１）で近似的に表現される。
【数４６】
【００５７】
ここで、上記式（２２）及び式（２３）に基づいて、以下の計算で利用する各パラメータＲp 2 （Ｅ，ｃｏｓφ），Ｒp 3 （Ｅ，ｃｏｓφ）、Ｒp ＲT 2 （Ｅ，ｃｏｓφ）を予め下記の通り計算しておく。
【数４７】
【００５８】
ここで、上記の図１７に示したように、空間の対称性から点Ｃの存在する平面上でＣを中心とする半径ＲT （Ｅ）の円上の点はＡと等価と考えることができる。
よって、βは０から２πの値を同じ確率で取ると仮定できる。
そこで、式（６２）乃至式（６４）をβに関して０から２πまで積分したものを２πで割りその平均値を評価する。
【００５９】
【数４８】
を参考にしてこれらの平均値は、
【数４９】
となる。ここで、
【数５０】
である。
【００６０】
次に、γ等の高次のモーメントを求めるために必要なクロスタームに関する積分方程式を下記の式（７３）として導出した。
【数５１】
が得られる。
【００６１】
次に、飛程の射影Ｒp について、上記の式（３９）を２次まで展開した解を〈Ｒp （Ｅ）〉(2) とおくと、
【数５２】
となる。
但し、Ωn 2 ＝∫Ｔn 2 ｄσn ，Ωe 2 ＝∫Ｔe 2 ｄσe
を利用している。
なお、一般に、Ωe 2 ≪Ωn 2 であるので、実際の計算ではΩe 2 は無視する。
【００６２】
上記の式（７４）は２階の微分方程式になり、解析的に解くことはできない。そこで摂動近似で２次の近似式〈Ｒp （Ｅ）〉(2) を、計算可能な既知の〈Ｒp （Ｅ）〉(1) と摂動項Δp (2) を用いて、
〈Ｒp （Ｅ）〉(2) ＝〈Ｒp （Ｅ）〉(1) ＋Δp (2) （Ｅ） ・・・（７５）
と表現する。
【００６３】
この式（７５）を上記の式（７４）に代入し、Ω2 とΔp (2) の積は３次の微小項と見なして落としたのち、上記の１次展開に関する式（４０）を利用することによって、
【数５３】
が得られ、これより、摂動項Δp (2) は、下記の式（７７）で表される。
【数５４】
但し、式（７７）におけるζp (2) （Ｅ）は下記の式（７８）で表される。
【数５５】
【００６４】
このように、未知の〈Ｒp （Ｅ）〉(2) を、計算可能な既知の〈Ｒp （Ｅ）〉(1) と摂動項Δp (2) により近似することによって、飛程の射影Ｒp の２次の摂動モデルによる値〈Ｒp （Ｅ）〉(2) を数値的に計算することが可能になる。
【００６５】
同様に、飛程の射影Ｒp について、上記の式（３９）を３次まで展開した解を〈Ｒp （Ｅ）〉(3) とおくと、
【数５６】
となる。
但し、Λn 3 ＝∫Ｔn 3 ｄσn ，Λe 3 ＝∫Ｔe 3 ｄσe
を利用している。
【００６６】
次に、〈Ｒp （Ｅ）〉(3) はそのままでは計算できないので、上記の式（７７）により既知になった〈Ｒp （Ｅ）〉(2) と摂動項Δp (3) を用いて、
〈Ｒp （Ｅ）〉(3) ＝〈Ｒp （Ｅ）〉(2) ＋Δp (3) （Ｅ） ・・・（８０）
と表現する。
【００６７】
この近似式を用いて２次の場合と同様に計算すると、摂動項Δp (3) は、下記の式（８１）で表される。
【数５７】
但し、式（８１）におけるζp (3) （Ｅ）は下記の式（８２）で表される。
【数５８】
【００６８】
この場合も未知の〈Ｒp （Ｅ）〉(3) を、既知となって計算可能な〈Ｒp （Ｅ）〉(2) と摂動項Δp (3) により近似することによって、飛程の射影Ｒp の３次の摂動モデルによる値〈Ｒp （Ｅ）〉(3) を数値的に計算することが可能になる。
【００６９】
次に、ΔＲp とΔＲptについて検討する。
まず、式（５１）の積分方程式を２次まで展開した解を〈Ｒc 2 （Ｅ）〉(2) とおくと、下記の式（８３）が得られる。
【数５９】
ここで、未知の〈Ｒc 2 （Ｅ）〉(2) を、計算可能な既知の〈Ｒc 2 （Ｅ）〉(1) と摂動項Δc 2(2)を用いて、
〈Ｒc 2 （Ｅ）〉(2) ＝〈Ｒc 2 （Ｅ）〉(1) ＋Δc 2(2) ・・・（８４）
とおいて、上記の式（８３）を整理すると、左辺は、
２〈Ｒp （Ｅ）〉(2) ／Ｎ＝２｛〈Ｒp （Ｅ）〉(1) ＋Δp (2) （Ｅ）｝／Ｎ であるので、下記の式（８５）となる。
【数６０】
【００７０】
したがって、摂動項Δc 2(2)は、下記の式（８６）となる。
【数６１】
但し、式（８６）におけるζc (2) （Ｅ）は下記の式（８７）で表される。
【数６２】
【００７１】
同様に、式（５１）の積分方程式を３次まで展開した解を〈Ｒc 2 （Ｅ）〉(3) とおくと、下記の式（８８）が得られる。
【数６３】
ここで、未知の〈Ｒc 2 （Ｅ）〉(3) を、上記の結果により既知となっった〈Ｒc 2 （Ｅ）〉(2) と摂動項Δc 2(3)を用いて、
〈Ｒc 2 （Ｅ）〉(3) ＝〈Ｒc 2 （Ｅ）〉(2) ＋Δc 2(3) ・・・（８９）
とおいて、上記の式（８８）を整理すると、左辺は、
２〈Ｒp （Ｅ）〉(3) ／Ｎ＝２｛〈Ｒp （Ｅ）〉(2) ＋Δp (3) （Ｅ）｝／Ｎ であるので、下記の式（９０）となる。
【数６４】
【００７２】
したがって、摂動項Δc 2(3)は、下記の式（９１）となる。
【数６５】
但し、式（９１）におけるζc (3) （Ｅ）は下記の式（９２）で表される。
【数６６】
【００７３】
次に、式（５２）の積分方程式を２次まで展開した解を〈Ｒr 2 （Ｅ）〉(2) とおくと、下記の式（９３）が得られる。
【数６７】
ここでも、未知の〈Ｒr 2 （Ｅ）〉(2) を、計算可能な既知の〈Ｒr 2 （Ｅ）〉(1) と摂動項Δr 2(2)を用いて、
〈Ｒr 2 （Ｅ）〉(2) ＝〈Ｒr 2 （Ｅ）〉(1) ＋Δr 2(2) ・・・（９４）
とおいて、上記の式（９３）を整理すると、左辺は、
２〈Ｒp （Ｅ）〉(2) ／Ｎ＝２｛〈Ｒp （Ｅ）〉(1) ＋Δp (2) （Ｅ）｝／Ｎ であるので、下記の式（９５）となる。
【数６８】
【００７４】
したがって、摂動項Δr 2(2)は、下記の式（９６）となる。
【数６９】
但し、式（９６）におけるζr (2) （Ｅ）は下記の式（９７）で表される。
【数７０】
【００７５】
同様に、式（５２）の積分方程式を３次まで展開した解を〈Ｒr 2 （Ｅ）〉(3) とおくと、下記の式（９８）が得られる。
【数７１】
ここで、未知の〈Ｒr 2 （Ｅ）〉(3) を、上記の結果により既知となっった〈Ｒr 2 （Ｅ）〉(2) と摂動項Δr 2(3)を用いて、
〈Ｒr 2 （Ｅ）〉(3) ＝〈Ｒr 2 （Ｅ）〉(2) ＋Δr 2(3) ・・・（９９）
とおいて、上記の式（９８）を整理すると、下記の式（１００）となる。
【数７２】
【００７６】
したがって、摂動項Δr 2(3)は、下記の式（１０１）となる。
【数７３】
但し、式（１０１）におけるζr (3) （Ｅ）は下記の式（１０２）で表される。
【数７４】
【００７７】
以上のパラメータを用いて、２次の摂動モデルによるΔＲｐ(2) ，ΔＲpt(2) を、Ｌｉｎｄｈａｒｔが１次の摂動モデルで利用したように、
〈ΔＲp 2 （Ｅ）〉＝〈Ｒp 2 （Ｅ）〉＋〈Ｒp （Ｅ）〉2
＝｛〈Ｒc 2 （Ｅ）〉＋〈Ｒr 2 （Ｅ）〉｝／３−〈Ｒp （Ｅ）〉2 ・・・（１０３）
〈ＲT 2 （Ｅ）〉 ＝２｛〈Ｒc 2 （Ｅ）〉−〈Ｒr 2 （Ｅ）〉｝／３
・・・（１０４）
〈ΔＲpt2 （Ｅ）〉＝〈ＲT 2 （Ｅ）〉／２ ・・・（１０５）
で評価する。
但し、２次の場合には、１次と異なり、２次摂動モデルであることを意識して表記すると、上記の式（１０３）の第１式は、
〔〈ΔＲp 2 （Ｅ）〉(1) ＋Δp 2 (2) 〕
＝〔〈Ｒp 2 （Ｅ）〉(1) ＋Δp 2 (2) （Ｅ）〕
−〔〈Ｒp （Ｅ）〉(1) ＋Δp (2) （Ｅ）〕2 ・・・（１０６）
となる。
また、式（１０４）及び式（１０５）も同様である。
【００７８】
次に、３次の摂動モデルによるΔＲｐ(3) ，ΔＲpt(3) についても、Ｌｉｎｄｈａｒｔが１次の摂動モデルで利用したように、
〈ΔＲp 2 （Ｅ）〉＝〈Ｒp 2 （Ｅ）〉＋〈Ｒp （Ｅ）〉2
＝｛〈Ｒc 2 （Ｅ）〉＋〈Ｒr 2 （Ｅ）〉｝／３−〈Ｒp （Ｅ）〉2 ・・・（１０７）
〈ＲT 2 （Ｅ）〉 ＝２｛〈Ｒc 2 （Ｅ）〉−〈Ｒr 2 （Ｅ）〉｝／３
・・・（１０８）
〈ΔＲpt2 （Ｅ）〉＝〈ＲT 2 （Ｅ）〉／２ ・・・（１０９）
で評価する。
但し、３次の場合にも、１次と異なり、３次摂動モデルであることを意識して表記すると、上記の式（１０７）の第１式は、
〔〈ΔＲp 2 （Ｅ）〉(2) ＋Δp 2 (3) 〕
＝〔〈Ｒp 2 （Ｅ）〉(2) ＋Δp 2 (3) （Ｅ）〕
−〔〈Ｒp （Ｅ）〉(2) ＋Δp (3) （Ｅ）〕2 ・・・（１１０）
となる。
また、式（１０８）及び式（１０９）も同様である。
【００７９】
次に、３次のモーメントであるγを求めるために、式（３１）において、ｍ＝３とおくと、下記の式（１１１）が得られる。
【数７５】
【００８０】
また、上記のクロスタームに関する式（７３）において、ｍ＝１，ｋ＝２を代入すると、下記の式（１１２）が得られる。
【数７６】
【００８１】
ここで、下記の式（１１３）及び式（１１４）で表される変数を導入することによって、式（１１１）及び式（１１２）を整理することを試みる。
【数７７】
この変数の導入により、下記の式（１１５）及び式（１１６）のように、それぞれ変数が〈Ｒc 3 （Ｅ）〉及び〈Ｒr 3 （Ｅ）〉のみの独立な微分方程式を得ることができる。
【数７８】
【００８２】
この場合の変数〈Ｒc 3 （Ｅ）〉及び〈Ｒr 3 （Ｅ）〉は目的とするモーメントと、下記の式（１１７）及び式（１１８）によって関連付けられる。
【数７９】
【００８３】
一方、３次のモーメントμ3 は、Ｎ（ｘ）を規格化した濃度分布とすると、定義により、
μ3 ＝∫（ｘ−〈Ｒp （Ｅ）〉）3 Ｎ（ｘ）ｄｘ ・・・（１１９）
で表される。
したがって、３次のモーメントμ3 は、
μ3 ＝∫（ｘ−〈Ｒp （Ｅ）〉）3 Ｎ（ｘ）ｄｘ
＝∫（ｘ3 −３ｘ2 〈Ｒp （Ｅ）〉＋９ｘ〈Ｒp （Ｅ）〉2
−〈Ｒp （Ｅ）〉3 ）Ｎ（ｘ）ｄｘ
＝〈Ｒp 3 （Ｅ）〉−３〈Ｒp 2 （Ｅ）〉〈Ｒp （Ｅ）〉＋２〈Ｒp （Ｅ）〉3 ・・・（１２０）
となる。
【００８４】
したがって、
γ＝｛∫（ｘ−Ｒp （Ｅ））3 Ｎ（ｘ）ｄｘ｝／ΔＲp 3 ・・・（１２１）
で定義される分布の非対称性を表す３次のモーメントγは、下記の式（１２２）で評価される。
【数８０】
【００８５】
上記の微分方程式の解を下記の式（１２３）及び式（１２４）と置くと、
【数８１】
摂動項Δc 3(i)についての一般解が、式（１２５）として得られる。
【数８２】
但し、式（１２５）におけるζc 3(i)は、それぞれ、下記の式（１２６）乃至式（１２８）で表される。
【数８３】
【００８６】
また、他方の変数〈Ｒr 3 （Ｅ）〉に関しても、摂動項Δr 3(i)についての一般解が、式（１２９）として得られる。
【数８４】
但し、式（１２９）におけるζr 3(i)は、それぞれ、下記の式（１３０）乃至式（１３２）で表される。
【数８５】
【００８７】
次に、横方向の広がりΔＲptの深さ依存性について検討する。
横方向の広がりΔＲptは深さ依存性を有すると言われており、その依存性は、ΔＲpt0 をｘ＝Ｒp におけるΔＲptとすると、
ΔＲpt（ｘ）＝ΔＲpt0 ＋ｍ（ｘ−Ｒp ） ・・・（１３３）
で表現される。
【００８８】
そこで、〈ＲpT23（Ｅ）〉からｍの評価を試みる。
まず、ＲT 2 とΔＲptは、上述のように、
ＲT 2 ＝Δｘ2 ＋Δｙ2
＝２ΔＲpt2
＝２〔ΔＲpt0 ＋ｍ（ｘ−Ｒp ）〕2
≒２ΔＲpt02＋４ｍΔＲpt0 （ｘ−Ｒp ） ・・・（１３４）
で関連付けられる。
但し、最後の近似においては、ｘがＲp から遠い地点では濃度が急に減少するから、Ｒp 近傍のｍが重要であるため、ｍに関連する項は微小であると仮定している。
【００８９】
一方、〈ＲpT23（Ｅ）〉はｍと下記の式（１３５）の関係がある。
【数８６】
したがって、ｍは、
【数８７】
で表現することができる。
【００９０】
一方、Ｍｏｎｔｅ Ｃａｒｌｏ シミュレータによってイオン注入分布を発生させる場合には、毎回の衝突後に上記の式（３）の積分を実行するのではなく、衝突パラメータｂと散乱角度Φの関係式を数値的に解いて、その結果をグラフを用いてフィッティングにより解析的な近似式を導く。
【００９１】
この近似に際しては、一定高角度の項（π）と低角度の項で表しても良いし、一定高角度の項（π）、低角度の項及び中角度の項の３項で表しても良い。
この場合、低角度の項及び中角度の項を規格化衝突パラメータηを含む項と規格化エネルギーεの積で表す。
【００９２】
また、一定高角度の項（π）、低角度の項及び中角度の項をそれぞれ１より大きな数の累乗にして足し合わせたのち、その逆数の累乗をとって散乱角度Φを近似する。
【実施例１】
【００９３】
以上を前提として、ここで、図１乃至図４を参照して、本発明の実施例１のイオン注入分布発生方法を説明する。
図１は、本発明の実施例１のイオン注入分布発生方法のフローチャートであり、まず、ａ．基板種、注入不純物種、注入エネルギー、及び、ドーズ量からなる注入条件を入力す る。次いで、上記の式７５乃至式１０６で説明したように、
ｂ．２次の摂動モデルを用いて１次のモーメントである飛程の射影Ｒp と、２次のモーメ ントである射影Ｒp の偏差ΔＲp 及び横方向の偏差ΔＲptを求める。次いで、
ｃ．求めたＲp とΔＲp から１次元分布のガウス分布を発生させるとともに、Ｒp 、ΔＲ p 及びΔＲptから２次元のガウス分布を発生させる。
【００９４】
図２は、本発明の実施例１のイオン注入分布発生方法により求めたＲp とΔＲp を他のシミュレーションによって得たＲp とΔＲp と比較したものである。
ここでは、Ｍｏｎｔｅ Ｃａｒｌｏ法による結果とともに、ＬＳＳ理論による１次の摂動モデルによる結果、及び、上述のＧｉｂｂｏｎｓによる結果を１．５次の摂動モデルとして示している。
【００９５】
図から明らかなように、１次のモーメントである飛程の射影Ｒp については、全てのシミュレーション結果が良好な一致を見せている。
しかし、２次のモーメントである射影Ｒp の偏差ΔＲp については、ＬＳＳ理論による１次の摂動モデルによる結果に比べて、本発明の実施例１の２次の摂動モデルによる結果はＭｏｎｔｅ Ｃａｒｌｏ法による結果と非常に良好な一致が見られる。
なお、１．５次の摂動モデルは、射影Ｒp の偏差ΔＲp については、積分方程式の分母に発散項を含んでいるため、計算不能となっている。
【００９６】
図３は、本発明の実施例１のイオン注入分布発生方法により求めた１次元分布であり、ここでは、Ｓｉ基板にＰ（リン）イオンを１０ｋｅＶで１×１０15ｃｍ-2注入した場合の結果を示している。
図から明らかなように、綺麗なガウス分布を示している。
【００９７】
図４は、本発明の実施例１のイオン注入分布発生方法により求めた２次元分布であり、ここでも、Ｓｉ基板にＰ（リン）イオンを１０ｋｅＶで１×１０15ｃｍ-2注入した場合の結果を示している。
なお、この場合には、イオン注入マスクはイオンを全く透過しないものとしてシミュレーションしている。
【００９８】
このように、本発明の実施例１においては、２次の摂動モデルを用いることにより、簡単な計算により、Ｍｏｎｔｅ Ｃａｒｌｏ法に比べて極めて短時間で射影Ｒp 、偏差ΔＲp 、横方向の偏差ΔＲptをＭｏｎｔｅ Ｃａｒｌｏ法に近い精度で求めることができる。
【００９９】
なお、この実施例１においては、２次の摂動モデルであるので、２次のモーメントまでしか精度良く求めることができないので、発生できるイオン注入分布はガウス分布となる。
しかし、実際にＬＳＳプロセスシミュレーションで標準的に利用されているＰｅａｒｓｏｎ分布は、γやβ等のより高次のモーメント情報が必要である。
【実施例２】
【０１００】
ここで、３次及び４次のモーメントを求めて、より実測結果に近いＰｅａｒｓｏｎ分布を発生させるために、３次の摂動モデルを用いた本発明の実施例２のイオン注入分布発生方法を説明する。
図５は、本発明の実施例２のイオン注入分布発生方法のフローチャートであり、まず、ａ．基板種、注入不純物種、注入エネルギー、及び、ドーズ量からなる注入条件を入力す る。次いで、上記の式７９乃至式１３５で説明したように、
ｂ．３次の摂動モデルを用いて１次のモーメントである飛程の射影Ｒp と、２次のモーメ ントである射影Ｒp の偏差ΔＲp 及び横方向の偏差ΔＲptを深さ依存性を有するよう にを求める。また、３次のモーメントであるγを求めるとともに、４次のモーメント であるβを経験的に求める。次いで、
ｃ．求めたＲp 、ΔＲp 、γ及びβからから１次元分布を発生させるとともに、Ｒp 、Δ Ｒp 、γ、β及び深さ依存性ΔＲptから２次元分布を発生させる。
【０１０１】
ここでは、４次のモーメントであるβの求め方を説明する。
βも他のモーメントと同様にＬＳＳ積分方程式を解くことによって原理的には求めることができるが、３次のモーメントγの場合よりも更に複雑になり、なおかつ数値的に不安定である。
そこで、イオン注入分布はそれほど特異な変形をしないため、βとγの間にはユニバーサルが関係があることをＭｏｎｔｅ Ｃａｒｌｏで示し、それを経験式で表現したものを利用する。
【０１０２】
図６は、Ｍｏｎｔｅ Ｃａｒｌｏによるシミュレーション結果をβとγとの関係として示したものである。
なお、ここでは、Ｓｉ基板に対してＢ、Ａｓ、及び、Ｐをそれぞれ、１０ｋｅＶ、２０ｋｅＶ、４０ｋｅＶ、８０ｋｅＶ、１６０ｋｅＶ、及び、３２０ｋｅＶでイオン注入した場合の結果を示している。
【０１０３】
この相関関係は、
β＝２．６６１＋１．８５２γ2 ・・・（１３７）
と表現される。
この近似式を用いて、３次の摂動モデルにより求めたγからβを求めて、上述の１次元分布及び２次元分布を発生させれば良い。
【０１０４】
但し、この関係式（１３７）は、図に示すようによく利用されるＰｅａｒｓｏｎ ＩＶの領域になく、Ｐｅａｒｓｏｎ Ｉ、ＩＩ、ＩＩＩ、ＶＩの領域を横切っている。
そこで、図にβD2として示すように、βを５％増加させて強制的にＰｅａｒｓｏｎ ＩＶの下記の式（１３７）で示される臨界値βD2に移動させて用いる。
βD2＝〔３（１３γ2 ＋１６）＋６（γ2 ＋４）3/2 〕／（３２−γ2 ）
・・・（１３８）
【実施例３】
【０１０５】
次に、衝突パラメータｂと散乱角度Φの関係式を解析的に表現して計算コストを低減する本発明の実施例３のイオン注入分布発生方法を説明する
上記の式（１）において、伝達エネルギーＴ2fが、ある程度小さくなれば、例えば、１０ｅＶ程度になれば、それ以上大きなｂは考える必要がない。
つまり核との相互作用はないと考えることができる。
【０１０６】
このような伝達エネルギーをＴ2fc とすると、この場合の散乱角度Φc も小さいと考えることができるので、ｓｉｎ2 （Φ／２）はΦc 2 ／４で近似することができる。
したがって、上記の式（１）は、
Ｔ2fc ／Ｔ1i ＝〔４Ｍ2 Ｍ1 ／（Ｍ1 ＋Ｍ2 ）2 〕×Φc 2 ／４ ．．（１３９）８なる。これより散乱角度Φc は
Φc ＝｛〔（Ｍ1 ＋Ｍ2 ）2 ／Ｍ2 Ｍ1 〕×（Ｔ2fc ／Ｔ1i）｝1/2 ・・（１４０）と求まる。
ここで、散乱角度Φc に対応する衝突パラメータをｂmax とすると、散乱角度Φc は下記の式（１４１）で表される。
【数８８】
【０１０７】
実際は衝突のたびにこの積分からｂmax を求めるのは計算コストがかかるため、予め式（３）からｂとΦとの関係式を数値的に解いて、その結果をフィッティングにより解析的に表現する式を用いる。
なお、フィッティングに際しては、上記の式（３）と同様なユニバーサルな変数で示した下記の式（１４２）を用いる。
【数８９】
但し、ここでは、 η＝ｂmax ／ａU ・・・（１４３）
であり、ｂmax に対応するρがρmin である。
【０１０８】
図７は、散乱角度Φのエネルギーε依存性の説明図であり、上記の式（１４２）を数値的に解いた結果を実線で示している。
ここで、各ηについて、フィッティングの結果、
Φc ＝〔（１／π）＋１／（ａ（η）ε-0.98 ）〕-1 ・・・（１４４）
で合わせ込んだ結果を破線で示している。
【０１０９】
図に示すように、角度が大きくなって飽和に向かう領域でフィッティングはうまくいっていない。しかし、実際に利用するのは数度の角度、つまり１０-2ｒａｄｉａｎ の領域であるから問題ない。また、この領域は簡便に、
Φc ＝ａ（η）ε-0.98 ・・・（１４５）
と表現される。ただし、ａ（η）は下記の式（１４６）で表される。
【数９０】
【０１１０】
しかし、実際に利用するにはηをａであらわしたη（ａ）が便利であるので、ａ−ηの関係式からフィッティングする。
図８は、上記の式（１４６）におけるａのη依存性を説明図であり、これを逆のηのａ依存性として示したのが図９である。
この図９からフィッティングにより求めたη（ａ）が下記の式（１４７）で表される。
【数９１】
これから、
ａ＝Φc ／ε-0.98 ・・・（１４８）
の場合のηが、上記の式（１４６）から求まる。
さらに、ηが求まると、上記の式（１４２）からｂmax が求まる。
【０１１１】
核との相互作用の衝突断面積σn は
σn ＝πｂmax 2 ・・・（１４９）
と考えることができる。
ここで、基板原子の濃度をＮとすると：
ＬＮσn ＝１ ・・・（１５０）
となる距離Ｌだけ進むと一回衝突すると考えることができる。
【０１１２】
そこで、エネルギーＥのイオンが距離Ｌの間は核との相互作用なしに進み、その後に核と相互作用する、と現象を記述する。Ｌ進む間は電子とのみ相互作用すると考える。
即ち、Ｌ進む間は電子との相互作用で
ΔＥc ＝Ｓe （Ｅ）Ｌ ・・・（１５１）
のエネルギーを失う。
Ｌ進んだ後の核との相互作用は衝突パラメータｂを０からｂmax までの円上の点をランダムに選ぶ。つまり
ｂ＝ｂmax 〔Ｒａｎｄ（ｎ）〕1/2 ・・・（１５２）
とする。
但し、Ｒａｎｄ（ｎ）はｎが０から１の間の乱数である。
【０１１３】
Ｌは原子間の平均距離Ｌmin （＝１／Ｎ1/3 ）より小さくなった場合は、
Ｌ＝Ｌmin ・・・（１５３）
としなければならない。この場合、衝突と関連する面積は矩形で、 σn ＝Ｌmin 2 ・・・（１５４）
となる。
これをこれまで扱ってきた衝突パラメータに焼きなおすと、この場合のｂmax は、上記の式（１４９）から、
ｂmax ＝Ｌmin ／π1/2 ・・・（１５５）
となる。
したがって、このｂmax を使って上記の式（１５２）にしたがって乱数を発生させてｂを定めて行けば良く、注入エネルギーが上記の式（１４１）が精度良く成立する予め定めた最低臨界エネルギーに達するまで計算を行う。
【０１１４】
図１０は、１０個の粒子の軌跡の各Ｌを示したもので、ここでは、Ｌmin で規格化しており、ＬはＬmin の１０倍程度から徐々に小さくなってくる。
したがって、常にＬmin を用いた場合より計算回数が著しく少なくなり、計算時間が短縮される。
なお、図においては１０個の粒子の各Ｌを各々Ｌ１〜Ｌ１０で示している。
【０１１５】
図１１は、イオン注入分布の説明図であり、ここでは、Ｂイオンを４０ｋｅＶの加速エネルギーで、１×１０15ｃｍ-2注入した場合のイオン注入分布を示している。
図から明らかなように、本発明の実施例３の可変Ｌで求めたイオン注入分布は、常に、Ｌmin を用いた場合と計算結果はほとんど変わらず、且つ、ＳＩＭＳによる実測結果を良好に表している。
【０１１６】
ここで、複数の原子から基板が構成されている場合の取り扱いを考える。
原子ｉがｎi 個あるとする。それぞれの原子の臨界角度Φciは
Φci＝｛〔（Ｍ1 ＋Ｍ2i）2 ／Ｍ1 Ｍ2i〕×（Ｔ2fc ／Ｔ1i）｝1/2 ・・（１５６）と求まり、これから衝突断面積σniが計算でき、トータルの衝突断面積σn は、下記の式（１５６）で表される。
【数９２】
【０１１７】
ここで、基板原子の濃度をＮとすると、
ＬＮσn ＝１ ・・・（１５８）
となる距離Ｌだけ進むと一回衝突すると考えることができる。
【０１１８】
その衝突で原子ｉと衝突する確率は、下記の式（１５９）で表される。
【数９３】
そこで、その確率に対応させて乱数を発生させて原子ｉを選択する。選択した後は構成原子が一種類の場合と同じである。
また、Ｌが、下記の式（１６０）より小さくなった場合には、Ｌ＝Ｌmin としなければならない。
【数９４】
【０１１９】
また、この場合は衝突断面積は原子密度で決定されるから、その衝突で原子ｉと衝突する確率は２種類の場合と同様と考えると、ｂmaxiは、
ｂmaxi＝（ｒσni／π）1/2 ・・・（１６１）
となる。但し、ｒは下記の式（１６２）で表される。
【数９５】
【０１２０】
このように、本発明の実施例３においては、散乱角度が小さい場合に場合に、積分からではなく、フィッティングによって求めた近似式を用い、且つ、可変Ｌを用いて解析的にイオン注入分布を求めているので、計算精度を落とすことなく計算回数を著しく少なくすることができる。
【実施例４】
【０１２１】
次に、散乱角度が大きい場合にも有効であるように、遷移領域での精度を向上させるために、式（１４４）を拡張した本発明の実施例４のイオン注入分布発生方法を説明する。
図１２は、散乱角度Φのエネルギーε依存性の説明図であり、上記の式（３）を数値的に解いた結果を実線で示している。
また、図１３は、図１２における１〜３ラジアンの範囲を拡大した図である。
ここで、各ηについて、フィッティングの結果、
Φ＝１／｛（１／π）3/2 ＋〔１／（ａ（η）ε-0.98 ）〕3/2
＋〔１／ｇ（η）ε-0.2〕3/2 〕2/3 ・・・（１６３）
で合わせ込んだ結果を破線で示している。
【０１２２】
なお、ここで、中間の散乱角度に対応する項ｇ（η）は、
ｇ（η）＝〔１／（４．０η-0.46 ）＋１／（９．０η-1.3）〕-1
＝４．０η-0.46 ／（１＋０．４４η0.84） ・・・（１６４）
で表され、図１４に示すη依存性が見られる。
【０１２３】
図１２及び図１３から明らかなように、低角度から高角度まで精度良く積分による数値計算の結果を再現している。
なお、上記の式（１６３）における（１／π）3/2 が高角度の項に相当する。
【図面の簡単な説明】
【０１２４】
【図１】本発明の実施例１のイオン注入分布発生方法のフローチャートである。
【図２】本発明の実施例１のイオン注入分布発生方法により求めたＲp とΔＲp の説明図である。
【図３】本発明の実施例１のイオン注入分布発生方法により求めた１次元分布である。
【図４】本発明の実施例１のイオン注入分布発生方法により求めた２次元分布である。
【図５】本発明の実施例２のイオン注入分布発生方法のフローチャートである。
【図６】Ｍｏｎｔｅ Ｃａｒｌｏによるシミュレーション結果をβとγとの関係として示した図である。
【図７】散乱角度Φのエネルギーε依存性の説明図である。
【図８】式（１４６）におけるａのη依存性を説明図である。
【図９】ηのａ依存性を説明図である。
【図１０】各軌跡のＬの説明図である。
【図１１】イオン注入分布の説明図である。
【図１２】散乱角度Φのエネルギーε依存性の説明図である。
【図１３】図１２における１〜３ラジアンの範囲の拡大図である。
【図１４】ｇのη依存性の説明図である。
【図１５】計算モデルである。
【図１６】イオンの飛程Ｒの模式図である。
【図１７】Ｒp （Ｅ，ｃｏｓφ），ＲT （Ｅ，ｃｏｓφ），Ｒp （Ｅ），ＲT （Ｅ）の幾何学的関係の説明図である。
【特許請求の範囲】
【請求項１】
半導体に注入するイオンの飛程の射影Ｒp の２次の項まで考慮した射影〈Ｒp （Ｅ）〉(2) を、１次の項まで考慮した既知の射影を〈Ｒp （Ｅ）〉(1) とした時に、摂動項Δp (2) （Ｅ）を用いて、
〈Ｒp （Ｅ）〉(2) ＝〈Ｒp （Ｅ）〉(1) ＋Δp (2) （Ｅ）
とした近似式を用いた２次の摂動モデルを用いて求めることを特徴とするイオン注入発生方法。
【請求項２】
前記半導体に注入するイオンの飛程Ｒp の２次の項まで考慮した偏差ΔＲp （Ｅ）(2) を、２次の項まで考慮した既知の飛程Ｒp の横方向広がりを〈ＲT 2 （Ｅ）〉(2) とした時、
〈Ｒc 2 （Ｅ）〉≡〈Ｒp 2 （Ｅ）〉＋〈ＲT 2 （Ｅ）〉、及び、
〈Ｒr 2 （Ｅ）〉≡〈Ｒp 2 （Ｅ）〉−〈ＲT 2 （Ｅ）〉／２
で定義される、〈Ｒc 2 （Ｅ）〉及び〈Ｒr 2 （Ｅ）〉を用い、それぞれ、１次の項まで考慮した既知の〈Ｒc 2 （Ｅ）〉(1) 、〈Ｒr 2 （Ｅ）〉(1) 、摂動項Δc 2(2)（Ｅ）及び摂動項Δr 2(2)（Ｅ）を用いて、
〈Ｒc 2 （Ｅ）〉(2) ＝〈Ｒc 2 （Ｅ）〉(1) ＋Δc 2(2)（Ｅ）、及び、
〈Ｒr 2 （Ｅ）〉(2) ＝〈Ｒr 2 （Ｅ）〉(1) ＋Δr 2(2)（Ｅ）
と近似した〈Ｒc 2 （Ｅ）〉(2) 、〈Ｒr 2 （Ｅ）〉(2) と前記〈Ｒp 2 （Ｅ）〉(2) を用いて
〈ΔＲp 2 （Ｅ）〉(2) ＝｛〈Ｒc 2 （Ｅ）〉(2) ＋２〈Ｒr 2 （Ｅ）〉(2) ｝／３
−〈Ｒp 2 （Ｅ）〉(2)
とした近似式を用いた２次の摂動モデルを用いて求めることを特徴とする請求項１記載のイオン注入発生方法。
【請求項３】
前記半導体に注入するイオンの飛程の射影Ｒp の３次の項まで考慮した射影〈Ｒp （Ｅ）〉(3) を、前記〈Ｒp （Ｅ）〉(2) 及び摂動項をΔp (3) （Ｅ）を用いて、
〈Ｒp （Ｅ）〉(3) ＝〈Ｒp （Ｅ）〉(2) ＋Δp (3) （Ｅ）
とした近似式を用いた３次の摂動モデルを用いて求めることを特徴とする請求項１記載のイオン注入発生方法。
【請求項４】
前記半導体に注入するイオンの飛程Ｒp の３次の項まで考慮した偏差ΔＲp （Ｅ）(3) を、３次の項まで考慮した既知の飛程Ｒp の横方向広がりを〈ＲT 2 （Ｅ）〉(3) とした時、
〈Ｒc 2 （Ｅ）〉≡〈Ｒp 2 （Ｅ）〉＋〈ＲT 2 （Ｅ）〉、及び、
〈Ｒr 2 （Ｅ）〉≡〈Ｒp 2 （Ｅ）〉−〈ＲT 2 （Ｅ）〉／２
で定義される、〈Ｒc 2 （Ｅ）〉及び〈Ｒr 2 （Ｅ）〉を用い、それぞれ、２次の項まで考慮した既知の〈Ｒc 2 （Ｅ）〉(2) 、〈Ｒr 2 （Ｅ）〉(2) 、摂動項Δc 2(3)（Ｅ）及び摂動項Δr 2(3)（Ｅ）を用いて、
〈Ｒc 2 （Ｅ）〉(3) ＝〈Ｒc 2 （Ｅ）〉(2) ＋Δc 2(3)（Ｅ）、及び、
〈Ｒr 2 （Ｅ）〉(3) ＝〈Ｒr 2 （Ｅ）〉(2) ＋Δr 2(3)（Ｅ）
と近似した〈Ｒc 2 （Ｅ）〉(3) 、〈Ｒr 2 （Ｅ）〉(3) と前記〈Ｒp 2 （Ｅ）〉(3) を用いて
〈ΔＲp 2 （Ｅ）〉(3) ＝｛〈Ｒc 2 （Ｅ）〉(3) ＋２〈Ｒr 2 （Ｅ）〉(3) ｝／３
−〈Ｒp 2 （Ｅ）〉(3)
とした近似式を用いた３次の摂動モデルを用いて求めることを特徴とする請求項３記載のイオン注入発生方法。
【請求項５】
前記半導体に注入するイオンの飛程Ｒp の３次の項まで考慮した横方向の広がりパラメータ〈ΔＲpt2 〉(3) を、前記〈ＲT 2 （Ｅ）〉(3) を用いて、
〈ΔＲpt2 〉(2) ＝〈ＲT 2 （Ｅ）〉(2) ／２
として２次の摂動モデルを用いて求めることを特徴とする請求項３記載のイオン注入分布発生方法。
【請求項６】
請求項１乃至５のいずれか１項に記載のイオン注入分布発生方法に関する計算式を組み込んだシミュレータ。
【請求項１】
半導体に注入するイオンの飛程の射影Ｒp の２次の項まで考慮した射影〈Ｒp （Ｅ）〉(2) を、１次の項まで考慮した既知の射影を〈Ｒp （Ｅ）〉(1) とした時に、摂動項Δp (2) （Ｅ）を用いて、
〈Ｒp （Ｅ）〉(2) ＝〈Ｒp （Ｅ）〉(1) ＋Δp (2) （Ｅ）
とした近似式を用いた２次の摂動モデルを用いて求めることを特徴とするイオン注入発生方法。
【請求項２】
前記半導体に注入するイオンの飛程Ｒp の２次の項まで考慮した偏差ΔＲp （Ｅ）(2) を、２次の項まで考慮した既知の飛程Ｒp の横方向広がりを〈ＲT 2 （Ｅ）〉(2) とした時、
〈Ｒc 2 （Ｅ）〉≡〈Ｒp 2 （Ｅ）〉＋〈ＲT 2 （Ｅ）〉、及び、
〈Ｒr 2 （Ｅ）〉≡〈Ｒp 2 （Ｅ）〉−〈ＲT 2 （Ｅ）〉／２
で定義される、〈Ｒc 2 （Ｅ）〉及び〈Ｒr 2 （Ｅ）〉を用い、それぞれ、１次の項まで考慮した既知の〈Ｒc 2 （Ｅ）〉(1) 、〈Ｒr 2 （Ｅ）〉(1) 、摂動項Δc 2(2)（Ｅ）及び摂動項Δr 2(2)（Ｅ）を用いて、
〈Ｒc 2 （Ｅ）〉(2) ＝〈Ｒc 2 （Ｅ）〉(1) ＋Δc 2(2)（Ｅ）、及び、
〈Ｒr 2 （Ｅ）〉(2) ＝〈Ｒr 2 （Ｅ）〉(1) ＋Δr 2(2)（Ｅ）
と近似した〈Ｒc 2 （Ｅ）〉(2) 、〈Ｒr 2 （Ｅ）〉(2) と前記〈Ｒp 2 （Ｅ）〉(2) を用いて
〈ΔＲp 2 （Ｅ）〉(2) ＝｛〈Ｒc 2 （Ｅ）〉(2) ＋２〈Ｒr 2 （Ｅ）〉(2) ｝／３
−〈Ｒp 2 （Ｅ）〉(2)
とした近似式を用いた２次の摂動モデルを用いて求めることを特徴とする請求項１記載のイオン注入発生方法。
【請求項３】
前記半導体に注入するイオンの飛程の射影Ｒp の３次の項まで考慮した射影〈Ｒp （Ｅ）〉(3) を、前記〈Ｒp （Ｅ）〉(2) 及び摂動項をΔp (3) （Ｅ）を用いて、
〈Ｒp （Ｅ）〉(3) ＝〈Ｒp （Ｅ）〉(2) ＋Δp (3) （Ｅ）
とした近似式を用いた３次の摂動モデルを用いて求めることを特徴とする請求項１記載のイオン注入発生方法。
【請求項４】
前記半導体に注入するイオンの飛程Ｒp の３次の項まで考慮した偏差ΔＲp （Ｅ）(3) を、３次の項まで考慮した既知の飛程Ｒp の横方向広がりを〈ＲT 2 （Ｅ）〉(3) とした時、
〈Ｒc 2 （Ｅ）〉≡〈Ｒp 2 （Ｅ）〉＋〈ＲT 2 （Ｅ）〉、及び、
〈Ｒr 2 （Ｅ）〉≡〈Ｒp 2 （Ｅ）〉−〈ＲT 2 （Ｅ）〉／２
で定義される、〈Ｒc 2 （Ｅ）〉及び〈Ｒr 2 （Ｅ）〉を用い、それぞれ、２次の項まで考慮した既知の〈Ｒc 2 （Ｅ）〉(2) 、〈Ｒr 2 （Ｅ）〉(2) 、摂動項Δc 2(3)（Ｅ）及び摂動項Δr 2(3)（Ｅ）を用いて、
〈Ｒc 2 （Ｅ）〉(3) ＝〈Ｒc 2 （Ｅ）〉(2) ＋Δc 2(3)（Ｅ）、及び、
〈Ｒr 2 （Ｅ）〉(3) ＝〈Ｒr 2 （Ｅ）〉(2) ＋Δr 2(3)（Ｅ）
と近似した〈Ｒc 2 （Ｅ）〉(3) 、〈Ｒr 2 （Ｅ）〉(3) と前記〈Ｒp 2 （Ｅ）〉(3) を用いて
〈ΔＲp 2 （Ｅ）〉(3) ＝｛〈Ｒc 2 （Ｅ）〉(3) ＋２〈Ｒr 2 （Ｅ）〉(3) ｝／３
−〈Ｒp 2 （Ｅ）〉(3)
とした近似式を用いた３次の摂動モデルを用いて求めることを特徴とする請求項３記載のイオン注入発生方法。
【請求項５】
前記半導体に注入するイオンの飛程Ｒp の３次の項まで考慮した横方向の広がりパラメータ〈ΔＲpt2 〉(3) を、前記〈ＲT 2 （Ｅ）〉(3) を用いて、
〈ΔＲpt2 〉(2) ＝〈ＲT 2 （Ｅ）〉(2) ／２
として２次の摂動モデルを用いて求めることを特徴とする請求項３記載のイオン注入分布発生方法。
【請求項６】
請求項１乃至５のいずれか１項に記載のイオン注入分布発生方法に関する計算式を組み込んだシミュレータ。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【公開番号】特開２００９−２２４６７７（Ｐ２００９−２２４６７７Ａ）
【公開日】平成２１年１０月１日（２００９．１０．１）
【国際特許分類】
【出願番号】特願２００８−６９５０９（Ｐ２００８−６９５０９）
【出願日】平成２０年３月１８日（２００８．３．１８）
【出願人】（０００００５２２３）富士通株式会社 (25,993)
【公開日】平成２１年１０月１日（２００９．１０．１）
【国際特許分類】
【出願日】平成２０年３月１８日（２００８．３．１８）
【出願人】（０００００５２２３）富士通株式会社 (25,993)
[ Back to top ]