(57)【要約】【課題】反応種の特性を考慮することで形状シミュレーションの高速化および高精度化を実現可能にする。【解決手段】一の実施形態による形状シミュレーション装置は、物質の表面を複数の計算要素に分割する分割部を備える。さらに、前記装置は、各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定する判定部を備える。さらに、前記装置は、前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する計算部を備える。さらに、前記判定部は、各計算要素に到達したイオン種が反射する場合における前記形態係数の計算用に前記判定を行う場合、前記イオン種の反射方向のカットオフ角度を設定し、前記直線を伸ばす方向を前記カットオフ角度の範囲内に制限して前記判定を行う。【選択図】図１７
2014135339.tif 000002
【特許請求の範囲】
【請求項１】
物質の表面を複数の計算要素に分割する分割部と、
各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定する判定部と、
前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する計算部とを備え、
前記判定部は、各計算要素に到達したイオン種が反射する場合における前記形態係数の計算用に前記判定を行う場合、前記イオン種の反射方向のカットオフ角度を設定し、前記直線を伸ばす方向を前記カットオフ角度の範囲内に制限して前記判定を行い、
前記判定部は、各計算要素に到達したイオン種によるスパッタにより中性種が発生する場合における前記形態係数の計算用に前記判定を行う場合、前記中性種の発生方向のカットオフ角度を設定し、前記直線を伸ばす方向を前記カットオフ角度の範囲内に制限して前記判定を行い、
前記判定部は、前記複数の計算要素のうちの第１の計算要素からの直線が第２の計算要素にぶつかった場合には、前記第１の計算要素からの直線が前記第２の計算要素の周囲の第３の計算要素にぶつかるか否かと、前記第３の計算要素が前記第１の計算要素の前記カットオフ角度の範囲内に位置するか否かとを判定し、
前記判定部は、前記第３の計算要素として、前記第２の計算要素に直接的に隣接する計算要素と、前記第２の計算要素に、肯定の判定結果が得られた計算要素を介して間接的に隣接する計算要素とを選択し、選択可能な前記第３の計算要素がなくなるまで前記判定を繰り返し、
前記判定部は、各計算要素に到達した中性種が各計算要素から再度飛散する場合における前記形態係数の計算用に前記判定を行う場合、前記直線を伸ばす方向のカットオフを行わずに前記判定を行う、
形状シミュレーション装置。
【請求項２】
物質の表面を複数の計算要素に分割する分割部と、
各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定する判定部と、
前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する計算部とを備え、
前記判定部は、各計算要素に到達したイオン種が反射する場合における前記形態係数の計算用に前記判定を行う場合、前記イオン種の反射方向のカットオフ角度を設定し、前記直線を伸ばす方向を前記カットオフ角度の範囲内に制限して前記判定を行い、
前記判定部は、前記複数の計算要素のうちの第１の計算要素からの直線が第２の計算要素にぶつかった場合には、前記第１の計算要素からの直線が前記第２の計算要素の周囲の第３の計算要素にぶつかるか否かと、前記第３の計算要素が前記第１の計算要素の前記カットオフ角度の範囲内に位置するか否かとを判定し、
前記判定部は、前記第３の計算要素として、前記第２の計算要素に直接的に隣接する計算要素と、前記第２の計算要素に、肯定の判定結果が得られた計算要素を介して間接的に隣接する計算要素とを選択し、選択可能な前記第３の計算要素がなくなるまで前記判定を繰り返す、
形状シミュレーション装置。
【請求項３】
前記判定部は、各計算要素に到達したイオン種によるスパッタにより中性種が発生する場合における前記形態係数の計算用に前記判定を行う場合、前記中性種の発生方向のカットオフ角度を設定し、前記直線を伸ばす方向を前記カットオフ角度の範囲内に制限して前記判定を行う、請求項２に記載の形状シミュレーション装置。
【請求項４】
前記判定部は、各計算要素に到達した中性種が各計算要素から再度飛散する場合における前記形態係数の計算用に前記判定を行う場合、前記直線を伸ばす方向のカットオフを行わずに前記判定を行う、請求項２または３に記載の形状シミュレーション装置。
【請求項５】
物質の表面を複数の計算要素に分割し、
各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定し、
前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する、
ことを含み、
前記判定では、各計算要素に到達したイオン種が反射する場合における前記形態係数の計算用に前記判定を行う場合、前記イオン種の反射方向のカットオフ角度を設定し、前記直線を伸ばす方向を前記カットオフ角度の範囲内に制限して前記判定を行い、
前記判定では、前記複数の計算要素のうちの第１の計算要素からの直線が第２の計算要素にぶつかった場合には、前記第１の計算要素からの直線が前記第２の計算要素の周囲の第３の計算要素にぶつかるか否かと、前記第３の計算要素が前記第１の計算要素の前記カットオフ角度の範囲内に位置するか否かとを判定し、
前記判定では、前記第３の計算要素として、前記第２の計算要素に直接的に隣接する計算要素と、前記第２の計算要素に、肯定の判定結果が得られた計算要素を介して間接的に隣接する計算要素とを選択し、選択可能な前記第３の計算要素がなくなるまで前記判定を繰り返す、
形状シミュレーション方法。
【請求項６】
形状シミュレーション方法をコンピュータに実行させる形状シミュレーションプログラムであって、
前記形状シミュレーション方法は、
物質の表面を複数の計算要素に分割し、
各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定し、
前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する、
ことを含み、
前記判定では、各計算要素に到達したイオン種が反射する場合における前記形態係数の計算用に前記判定を行う場合、前記イオン種の反射方向のカットオフ角度を設定し、前記直線を伸ばす方向を前記カットオフ角度の範囲内に制限して前記判定を行い、
前記判定では、前記複数の計算要素のうちの第１の計算要素からの直線が第２の計算要素にぶつかった場合には、前記第１の計算要素からの直線が前記第２の計算要素の周囲の第３の計算要素にぶつかるか否かと、前記第３の計算要素が前記第１の計算要素の前記カットオフ角度の範囲内に位置するか否かとを判定し、
前記判定では、前記第３の計算要素として、前記第２の計算要素に直接的に隣接する計算要素と、前記第２の計算要素に、肯定の判定結果が得られた計算要素を介して間接的に隣接する計算要素とを選択し、選択可能な前記第３の計算要素がなくなるまで前記判定を繰り返す、
形状シミュレーションプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明の実施形態は、形状シミュレーション装置、形状シミュレーション方法、および形状シミュレーションプログラムに関する。
【背景技術】
【０００２】
ＣＶＤ（Chemical Vapor Deposition）やＲＩＥ（Reactive Ion Etching）等による物質の表面の加工において、加工形状のシミュレーションは重要な技術となっている。このシミュレーションでは、物質の表面を複数の計算要素に分割し、各計算要素に到達する反応種のフラックスや、物質の局所的な表面成長速度を算出することが一般的である。しかしながら、フラックスや表面成長速度をすべての表面にて矛盾なく計算するには長い計算時間が必要となる。理由は、計算時間が計算要素の個数の２乗のオーダーで増えていくためである。また、反応種は、直進性が強く、異方的に入射するイオン種と、直進性が弱く、等方的に入射する中性種とに分類されるが、従来のシミュレーションはこれらの反応種の違いを考慮せずに行われているため、計算の無駄や誤差が生じている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２００９−４９１１０号公報
【非特許文献】
【０００４】
【非特許文献１】G. Kokkoris et al., "Simulation of SiO2 and Si feature etching for microelectronics and microelectromechanical systems fabrication: A combined simulator coupling modules of surface etching, local flux calculation, and profile evolution", J. Vac. Sci. Technol. A 22, 1896 (2004).
【非特許文献２】Arpan P. Mahorowala et al., "Etching of polysilicon in inductively coupled Cl2 and HBr discharges. II. Simulation of profile evolution using cellular representation of feature composition and Monte Carlo computation of flux and surface kinetics", J. Vac. Sci. Technol. B 20, 1064 (2002).
【非特許文献３】Shahram Abdollahi-Alibeik et al., "Analytical modeling of silicon etch process in high density plasma", J. Vac. Sci. Technol. A 17, 2485 (1999).
【発明の概要】
【発明が解決しようとする課題】
【０００５】
反応種の特性を考慮することで形状シミュレーションを高速化および高精度化することが可能な形状シミュレーション装置、形状シミュレーション方法、および形状シミュレーションプログラムを提供する。
【課題を解決するための手段】
【０００６】
一の実施形態による形状シミュレーション装置は、物質の表面を複数の計算要素に分割する分割部を備える。さらに、前記装置は、各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定する判定部を備える。さらに、前記装置は、前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する計算部を備える。さらに、前記判定部は、各計算要素に到達したイオン種が反射する場合における前記形態係数の計算用に前記判定を行う場合、前記イオン種の反射方向のカットオフ角度を設定し、前記直線を伸ばす方向を前記カットオフ角度の範囲内に制限して前記判定を行う。さらに、前記判定部は、前記複数の計算要素のうちの第１の計算要素からの直線が第２の計算要素にぶつかった場合には、前記第１の計算要素からの直線が前記第２の計算要素の周囲の第３の計算要素にぶつかるか否かと、前記第３の計算要素が前記第１の計算要素の前記カットオフ角度の範囲内に位置するか否かとを判定する。さらに、前記判定部は、前記第３の計算要素として、前記第２の計算要素に直接的に隣接する計算要素と、前記第２の計算要素に、肯定の判定結果が得られた計算要素を介して間接的に隣接する計算要素とを選択し、選択可能な前記第３の計算要素がなくなるまで前記判定を繰り返す。
【図面の簡単な説明】
【０００７】
【図１】第１実施形態の形状シミュレーション方法の手順を示したフローチャート図である。
【図２】第１実施形態における物質の初期構造の例を示す斜視図である。
【図３】レベルセット関数について説明するための模式図である。
【図４】図１のステップＳ３の詳細を示したフローチャート図である。
【図５】物質表面を複数の計算要素に分割した様子を示した模式図である。
【図６】イオン種と中性種の直進性の違いについて説明するための模式図である。
【図７】イオン種の反射方向のカットオフ角度について説明するための図である。
【図８】図４のステップＳ１２、Ｓ１３の詳細を示したフローチャート図である。
【図９】ローカル座標系について説明するための図である。
【図１０】可視判定値について説明するための模式図である。
【図１１】可視係数について説明するための模式図である。
【図１２】入射角度θinについて説明するための模式図である。
【図１３】鏡面境界条件について説明するための模式図である。
【図１４】周期境界条件について説明するための模式図である。
【図１５】２次元における計算要素可視判定値について説明するための模式図である。
【図１６】３次元における計算要素可視判定値について説明するための模式図である。
【図１７】図４のステップＳ１４の詳細を示したフローチャート図である。
【図１８】図１７の処理を説明するための模式図である。
【図１９】図１７の処理を説明するための模式図である。
【図２０】図１７の処理を説明するための模式図である。
【図２１】比較例における計算時間の例を示したグラフである。
【図２２】第１実施形態における計算時間の例を示したグラフである。
【図２３】第１実施形態と比較例の計算時間を比較したグラフである。
【図２４】第１実施形態と比較例におけるθ分割数と計算誤差との関係を示したグラフである。
【図２５】イオン種のみを取り扱った場合の第１実施形態と比較例の計算時間を比較したグラフである。
【図２６】イオン種と中性種を取り扱った場合の第１実施形態と比較例の計算時間を比較したグラフである。
【図２７】イオン種のみを取り扱った場合の第１実施形態と比較例の計算時間を項目別にに比較したグラフである。
【図２８】イオン種と中性種を取り扱った場合の第１実施形態と比較例の計算時間を項目別に比較したグラフである。
【図２９】第２実施形態の形状シミュレーション装置の構成を示す外観図である。
【図３０】図２９の制御部の構成を示すブロック図である。
【発明を実施するための最良の形態】
【０００８】
以下、本発明の実施形態を、図面を参照して説明する。これらの図面では、同一または類似の構成要素に同一の符号を付しており、重複する説明は適宜省略する。
【０００９】
（第１実施形態）
図１は、第１実施形態の形状シミュレーション方法の手順を示したフローチャート図である。本実施形態の形状シミュレーション方法は、パーソナルコンピュータやワークステーション等の情報処理装置を使用して行われる。
【００１０】
本実施形態の形状シミュレーション方法ではまず、物質の初期構造を情報処理装置に入力する（ステップＳ１）。図２は、第１実施形態における物質の初期構造の例を示す斜視図である。図２に示す初期構造は、シリコン基板１と、シリコン基板１上に順に形成されたシリコン窒化膜２およびシリコン酸化膜３と、シリコン窒化膜２とシリコン酸化膜３とを貫通する複数の貫通孔４を含んでいる。初期構造の入力方法の例としては様々なフォーマットが考えられるが、本実施形態では、物質の表面の形状を点列で表現し、情報処理装置がこれを読み取る方法を用いている。
【００１１】
次に、入力された初期構造から、初期レベルセット関数を作成する（ステップＳ２）。図３は、レベルセット関数について説明するための模式図である。レベルセット関数ψは、物質の表面からの距離ｄを用いて定義される関数であり、計算領域内のメッシュごとに値を有する。レベルセット関数ψの値は、物質の表面において０と定義される（ψ＝０）。また、物質の外部（真空中）ではψ＞０であり、物質の内部（物質中）ではψ＜０である。初期レベルセット関数を作成する際には、各メッシュ点から最近接となる表面を探し、その距離ｄを計算し、メッシュ点が真空中であればその符号を正とし、物質中であればその符号を負とする。なお、初期レベルセット関数は、ステップＳ２で作成する代わりに、ステップＳ１で入力してもよい。
【００１２】
次に、物質の局所的な表面成長速度Ｆを計算する（ステップＳ３）。ここで、表面の成長とは、表面への堆積だけでなく、表面のエッチングも含むものとする。なお、表面成長速度Ｆの計算は、タイムステップごとに行う必要はない。また、本実施形態では、後述するように、物質の表面におけるフラックス（総フラックス）から表面成長速度Ｆを計算し、表面成長速度Ｆからレベルセット関数を計算するが、代わりにフラックスからレベルセット関数を計算し、表面成長速度Ｆの計算は省略してもよい。
【００１３】
次に、表面成長速度Ｆを用いて、時間Δｔ経過後のレベルセット関数を計算する（ステップＳ４）。時間ｔにおけるレベルセット関数ψtは、以下の式（１）から計算できる。
【数１】
2014135339.tif 000003
【００１４】
ただし、∇はベクトル微分演算子を表し、｜∇ψt｜は∇ψtのノルムを表す。時間Δｔ経過後のレベルセット関数は、式（１）を離散化した式に従い、レベルセット関数を時間発展させることで計算可能である。なお、本実施形態では、表面形状を時間発展させる代わりに、ある表面形状における表面成長速度Ｆやフラックスを計算してもよい。これは、後述するステップＳ５を１ステップ目でＹｅｓと判定する場合に相当する。
【００１５】
次に、予め設定したプロセス時間が経過したか否かを判定する（ステップＳ５）。プロセス時間が終了した場合には、物質の最終形状を出力し（ステップＳ６）、計算終了となる。プロセス時間が終了していない場合には、ステップＳ３に戻る。
【００１６】
なお、本実施形態では、形状表現の手法としてレベルセット法を用いているが、レベルセット法以外のセル法やストリング法などの手法を用いてもよい。
【００１７】
（１）ステップＳ３の詳細
次に、図４を参照し、ステップＳ３の詳細について説明する。
【００１８】
図４は、図１のステップＳ３の詳細を示したフローチャート図である。
【００１９】
まず、レベルセット法で表された物質表面を、複数の計算要素に分割する（ステップＳ１１）。図５は、物質表面を複数の計算要素に分割した様子を示した模式図である。図５の例では、物質表面がメッシュごとに分割されている。その結果、１つのメッシュ内の物質表面が、１つの計算要素となっている。ステップＳ１１の処理を行うブロックは、本開示の分割部の例である。
【００２０】
なお、物質表面の分割方法は、メッシュ単位に限られるものではなく、どのような方法を採用してもよい。また、物質表面の分割は、タイムステップごとに行う必要はなく、例えばステップＳ１の直後に行ってもよい。
【００２１】
また、図５に示す計算領域は、２次元領域となっているが、代わりに３次元領域としてもよい。また、図５に示す各計算要素の形状は、線分となっているが、代わりに点や多角形などとしてもよい。
【００２２】
図５は、第１の計算要素ａと、第２の計算要素Ｂを示している。計算要素Ｂに到達する反応種のフラックスを計算する際には、気層から計算要素Ｂに直接的に到達する反応種のフラックスと、気層から任意の計算要素ａを介して計算要素Ｂに間接的に到達する反応種のフラックスの両方を考慮するのが一般的である。前者のフラックスを直接フラックスと呼び、後者のフラックスを間接フラックスと呼ぶ。また、これらの合計を、総フラックスと呼ぶ。なお、反応種の例としては、堆積種やエッチング種などが挙げられる。
【００２３】
また、反応種は、直進性が強く、異方的に入射するイオン種と、直進性が弱く、等方的に入射する中性種とに分類される。図６は、イオン種と中性種の直進性の違いについて説明するための模式図である。
【００２４】
図６は、計算要素ａに入射したイオン種が反射する様子を示している。図６に示すように、イオン種は直進性が強く、全方向に反射される訳ではないため、イオン種の間接フラックスを計算する際には、イオン種の反射方向のカットオフ角度を設定し、カットオフ角度の範囲外への反射を無視することが望ましい。これにより、本実施形態では、計算の無駄を減らして計算時間を短縮することや、無駄な計算の代わりに有用な計算に時間を割いて計算誤差を低減することができ、形状シミュレーションを高速化および高精度化することが可能となる。
【００２５】
図７は、イオン種の反射方向のカットオフ角度について説明するための図である。符号Ｅin、θinはそれぞれ、垂直入射したイオン種の入射方向と入射角度を示し、符号Ｅout、θoutはそれぞれ、垂直入射したイオン種を鏡面反射させた場合の反射方向と反射角度を示す。入射角度θinと反射角度θoutとの間には、θin＝θoutの関係が成り立つ。
【００２６】
本実施形態では、図７に示すように、反射方向Ｅoutを中心にカットオフ角度θcutを設定し、カットオフ角度θcutの範囲外へはイオン種は反射しないと想定する。なお、カットオフ角度θcutを使用した形状シミュレーションの詳細については、後述する。
【００２７】
以下、計算要素Ｂにおける中性種の総フラックスΓB,neとイオン種の総フラックスΓB,ionについて説明する。
【００２８】
計算要素Ｂにおける中性種の総フラックスΓB,neは、以下の式（２）のように、計算要素Ｂにおける中性種の直接フラックスΓB,ne-directと、任意の計算要素ａからの中性種の間接フラックスΓaB,ne-indirectの合計との和で表される。
【数２】
2014135339.tif 000004
【００２９】
同様に、計算要素Ｂにおけるイオン種の総フラックスΓB,ionは、以下の式（３）のように、計算要素Ｂにおけるイオン種の直接フラックスΓB,ion-directと、任意の計算要素ａからのイオン種の間接フラックスΓaB,ion-indirectの合計との和で表される。
【数３】
2014135339.tif 000005
【００３０】
ここで、間接フラックスΓaB,ne-indirect、ΓaB,ion-indirectは、例えば以下の式（４）、式（５）で表すことができる。
【数４】
2014135339.tif 000006
【数５】
2014135339.tif 000007
【００３１】
Ｓa(Γa,ion,Γa,ne)は、計算要素ａで吸収される中性種フラックスの割合を示す付着確率を表す。Ｒa(Γa,ion,Γa,ne)は、計算要素ａで反射されるイオン種フラックスの割合を示す反射確率を表す。Ｐa(Γa,ion,Γa,ne)は、計算要素ａでイオン種フラックスによるスパッタにより物質表面が削れて中性種フラックスが発生する割合を示すスパッタ確率を表す。Ｓa(Γa,ion,Γa,ne)、Ｒa(Γa,ion,Γa,ne)、Ｐa(Γa,ion,Γa,ne)の値は、計算要素ａにおけるイオン種総フラックスΓa,ionと中性種総フラックスΓa,neに依存する。
【００３２】
また、ν(ａ,Ｂ)は、計算要素ａと計算要素Ｂが互いに見えるか否かを示す可視係数（面間可視係数）を表す。計算要素ａ、Ｂを結ぶ直線が、計算要素ａ、Ｂ間において物質表面とぶつかる場合には、ν＝０となり、ぶつからない場合にはν＝１となる。
【００３３】
また、ｇ(ａ,Ｂ)は、計算要素ａと計算要素Ｂとの位置関係（面関係）を示す形態係数を表す。ｇ(ａ,Ｂ)の値は、計算要素ａ、Ｂの互いの見えやすさの程度を示している。ｇ(ａ,Ｂ)の値は、計算要素ａ、Ｂ間の距離や角度などに依存する。
【００３４】
形態係数ｇは、イオン種の反射やイオン種によるスパッタを取り扱う場合、イオン種の直進性や散乱にも依存する。そこで、本実施形態では、形態係数ｇの他に、イオン種の反射用の形態係数（反射形態係数）ｇionRと、イオン種によるスパッタ用の形態係数（スパッタ形態係数）ｇionSとを導入する。
【００３５】
ｇionR(ａ,Ｂ)は、計算要素ａに到達したイオン種が反射する場合における、計算要素ａと計算要素Ｂとの間の形態係数を表す。ｇionS(ａ,Ｂ)は、計算要素ａに到達したイオン種によるスパッタにより中性種が発生する場合における、計算要素ａと計算要素Ｂとの間の形態係数を表す。一方、ｇ(ａ,Ｂ)は、計算要素ａに到達した中性種が計算要素ａから再度飛散する場合等における、計算要素ａと計算要素Ｂとの間の形態係数を表している。中性種が再度飛散する場合の例としては、吸収された中性種が放出される場合や、中性種が反射する場合などが挙げられる。
【００３６】
式（２）、式（３）にそれぞれ式（４）、式（５）を代入すると、計算要素Ｂにおける総フラックスΓa,ion、Γa,neはそれぞれ、以下の式（６）、式（７）で表すことができる。
【数６】
2014135339.tif 000008
【数７】
2014135339.tif 000009
【００３７】
図４のフローでは次に、任意の計算要素における直接フラックスΓB,ne-direct、ΓB,ion-directや、任意の計算要素間における可視係数ν、形態係数ｇ、反射形態係数ｇionR、スパッタ形態係数ｇionSを計算する（ステップＳ１２〜Ｓ１４）。
【００３８】
次に、各計算要素ｉの直接フラックスΓi,ne-direct、Γi,ion-directをそれぞれ仮の総フラックスΓi,ne、Γi,ionとして使用し、各計算要素ｉにおける付着確率Ｓi(Γi,ion,Γi,ne)、反射確率Ｒi(Γi,ion,Γi,ne)、スパッタ確率Ｐi(Γi,ion,Γi,ne)を計算する（ステップＳ１５）。
【００３９】
次に、可視係数ν、形態係数ｇ、ｇionR、ｇionS、直接フラックスΓi,ne-direct、Γi,ion-direct、付着確率Ｓi(Γi,ion,Γi,ne)、反射確率Ｒi(Γi,ion,Γi,ne)、スパッタ確率Ｐi(Γi,ion,Γi,ne)を用いて、以下の式（８）、式（９）から、各計算要素ｉにおける総フラックスΓi,ion、Γi,neを計算する（ステップＳ１６）。
【数８】
2014135339.tif 000010
【数９】
2014135339.tif 000011
【００４０】
次に、ステップＳ１５とステップＳ１６の処理を、付着確率Ｓi(Γi,ion,Γi,ne)、反射確率Ｒi(Γi,ion,Γi,ne)、スパッタ確率Ｐi(Γi,ion,Γi,ne)の値が収束するまで繰り返す（ステップＳ１７）。なお、２回目以降のステップＳ１５では、前回のステップＳ１６で計算した総フラックスΓi,ion、Γi,neを、仮の総フラックスΓi,ion、Γi,neとして使用する。また、ステップＳ１７では、Ｓi(Γi,ion,Γi,ne)、Ｒi(Γi,ion,Γi,ne)、Ｐi(Γi,ion,Γi,ne)の値が収束したか否かを、Ｓi(Γi,ion,Γi,ne)、Ｒi(Γi,ion,Γi,ne)、Ｐi(Γi,ion,Γi,ne)の変化が閾値以下になったか否かで判定する。そして、これらの確率Ｓi(Γi,ion,Γi,ne)、Ｒi(Γi,ion,Γi,ne)、Ｐi(Γi,ion,Γi,ne)の値が収束した際の総フラックスΓi,ion、Γi,neを、総フラックスΓi,ion、Γi,neの正しい計算結果として取り扱う。
【００４１】
なお、計算要素の個数をＮとする場合、任意の計算要素間における可視係数νや形態係数ｇ、ｇionR、ｇionSは、まとめてＮ×Ｎ行列で表すことができる。行列の形で表された可視係数ν、形態係数ｇ、ｇionR、ｇionSを、それぞれ可視係数行列、形態係数行列と呼ぶ。また、任意の計算要素におけるフラックスは、まとめてＮ行ベクトルで表すことができる。ベクトルの形で表されたフラックスを、フラックスベクトルと呼ぶ。
【００４２】
この場合、式（８）は、以下の式（１０）のように行列方程式で表現することができる。
【数１０】
2014135339.tif 000012
【数１１】
2014135339.tif 000013
【数１２】
2014135339.tif 000014
【数１３】
2014135339.tif 000015
【００４３】
ただし、Ｉ、Ｊは処理対象の計算要素の個数を表し、例えばＩ＝Ｊ＝Ｎである。また、式（１３）では、表記スペースの都合上、付着確率Ｓj(Γj,ion,Γj,ne)をＳj(Γj)と略記すると共に、スパッタ確率Ｐi(Γi,ion,Γi,ne)を含む項の記載を省略している。
【００４４】
行列方程式（１０）は、どのような解法で解いてもよい。解法の例としては、反復法（ガウスザイデル法、ＳＯＲ法、ヤコビ法、共役勾配法など）や、直接法（ガウスの消去法、ＬＵ分解法、コレスキー分解法など）が挙げられる。行列方程式（１０）を解く際に、行列Ａneが疎行列である場合には、ＣＲＳなどの保存方法を用いた上で疎行列に適したルーチンを使用することで、計算処理の省メモリ化や高速化を図ってもよい。
【００４５】
なお、式（９）も、式（８）と同様に行列方程式で表現可能である。本実施形態では、これら２つの行列方程式を上記の解法で解くことができる。
【００４６】
図４のフローでは次に、総フラックスΓi,ion、Γi,neから、各計算要素ｉにおける局所的な表面成長速度Ｆiを計算する（ステップＳ１８）。例えば、Ｋ種類の反応種を使用する場合においては、表面成長速度Ｆiは、Ｋ個の局所的な総フラックスΓ1,i〜ΓK,iに依存した以下の式（１４）のような形でモデリングされる。
【数１４】
2014135339.tif 000016
【００４７】
ただし、ｋは１≦ｋ≦Ｋを満たす任意の実数である。なお、Ｋ種類の反応種は、中性種とイオン種のいずれか一方のみを含んでいてもよいし、中性種とイオン種の両方を含んでいてもよい。以上のようにして、ステップＳ３の処理が終了する。
【００４８】
（２）ステップＳ１２、Ｓ１３の詳細
次に、図８を参照し、ステップＳ１２、Ｓ１３の詳細について説明する。
【００４９】
ステップＳ１２、Ｓ１３では、中性種の直接フラックスΓB,ne-direct、イオン種の直接フラックスΓB,ion-direct、可視係数ν、および形態係数ｇを計算する。この際、中性種の直接フラックスΓB,ne-directと、イオン種の直接フラックスΓB,ion-directは、同じ方法で計算される。そこで、以下の説明では、中性種の直接フラックスΓB,ne-direct、可視係数ν、および形態係数ｇの計算方法について説明し、イオン種の直接フラックスΓB,ion-directの計算方法の説明は省略する。また、中性種の直接フラックスΓB,ne-directは、説明の便宜上、単に直接フラックスΓB,directと表記する。
【００５０】
図８は、図４のステップＳ１２、Ｓ１３の詳細を示したフローチャート図である。
【００５１】
図８のフローでは、各計算要素に固有のローカル座標系を使用する。図９は、ローカル座標系について説明するための図である。図９(ａ)は、各計算要素の法線ベクトルを示し、図９(ｂ)は、各計算要素におけるローカル座標系を示す。図９(ｂ)に示すように、ローカル座標系の直交座標（ｘlocal,ｙlocal,ｚlocal）は、＋ｚlocal方向が法線ベクトル方向と一致するように定められる。また、ローカル座標系の極座標（ｒlocal,θlocal,φlocal）は、極角θlocalが動径ｒlocalと＋ｚlocal方向との間の角度となり、偏角φlocalが動径ｒlocalと＋ｘlocal方向との間の角度となるように定められる。
【００５２】
計算要素Ｂにおける直接フラックスΓB,directは、以下の式（１５）により計算することができる。
【数１５】
2014135339.tif 000017
【００５３】
ただし、η(θlocal,φlocal)は、計算要素Ｂからθlocal,φlocalの方向に直線を伸ばした場合の可視判定の結果を示しており、可視判定値と呼ぶことにする。図１０は、可視判定値ηについて説明するための模式図である。図１０に示すように、上記直線が物質表面とぶつかる場合には、η＝０となり、ぶつからない場合にはη＝１となる。なお、図１０に示すように、物質表面の片側の方向にだけ直線を伸ばす場合には、式（１５）におけるθlocalの積分範囲は、０からπではなく、０からπ／２にしてもよい。
【００５４】
可視判定値ηと可視係数νとの違いについては、図１１を参照されたい。図１１は、可視係数νについて説明するための模式図である。ν(ａ,Ｂ)は、計算要素ａと計算要素Ｂが互いに見えるか否かを示す。計算要素ａ、Ｂを結ぶ直線が、計算要素ａ、Ｂ間において物質表面とぶつかる場合には、ν＝０となり、ぶつからない場合にはν＝１となる。前者の例としては計算要素ｄを、後者の例としては計算要素ｃを参照されたい。
【００５５】
また、ｆflatは、平坦面での直接フラックスを示し、入力値として事前に与えられる。さらに、Normは、以下の式（１６）で与えられる規格化定数を示す。さらに、ｆ(θlocal)は、直接フラックスの面積素片の係数を示し、例えば以下の式（１７）で与えられる。
【数１６】
2014135339.tif 000018
【数１７】
2014135339.tif 000019
【００５６】
ただし、θinは、図１２に示すような入射角度である。図１２は、入射角度θinについて説明するための模式図である。入射角度θinは、法線ベクトル方向とθlocal,φlocalの方向との間の角度に相当する。よって、ローカル座標系（ｒlocal,θlocal,φlocal）を用いる場合には、θin＝θlocalが成り立つ。
【００５７】
以下、図８のフローについて具体的に説明する。
【００５８】
まず、極角θlocalの数列θlocal(ｍ)の値（ｍ＝０、１、…、Ｍ−１）と、偏角φlocalの数列φlocal(ｏ)の値（ｏ＝０、１、…、Ｏ−１）とを計算する（ステップＳ２１）。これは、０からπまでの極角θlocalの範囲をＭ個の領域に分割し、０から２πまでの偏角φlocalの範囲をＯ個の領域に分割することに相当する。後述するように、式（１５）の積分計算は、これらの数列θlocal(ｍ),φlocal(ｏ)を用いて離散化される。
【００５９】
式（１５）の直接フラックスΓB,directの計算に、式（１７）に示す面積素片係数を使用する場合には、例えば、以下の式（１８）、式（１９）のような数列θlocal(ｍ),φlocal(ｏ)を用意する。
【数１８】
2014135339.tif 000020
【数１９】
2014135339.tif 000021
【００６０】
ここで、数列∂(ｍ)は、以下の式（２０）で与えられる。
【数２０】
2014135339.tif 000022
【００６１】
式（１８）のθlocal(ｍ)は、ｆ(θlocal)｜sinθlocal｜をθlocal＝０からθlocal＝θlocal(ｍ)まで積分した場合に、積分結果が∂(ｍ)となる角度を表す。この定義から式（２１）の関係が成り立ち、式（２１）から式（２２）が導出され、式（２２）を変形することで式（１８）が得られる。
【数２１】
2014135339.tif 000023
【数２２】
2014135339.tif 000024
【００６２】
以上のように、ステップＳ２１では、０からπまでの極角θlocalの範囲を非等間隔に分割し、０から２πまでの偏角φlocalの範囲を等間隔に分割する。なお、本実施形態では、極角θlocalの範囲だけでなく、偏角φlocalの範囲も非等間隔に分割してもよい。また、極角θlocalの積分範囲を０からπ／２とする場合には、０からπまでではなく、０からπ／２までの極角θlocalの範囲をＭ個の領域に分割するようにしてもよい。
【００６３】
次に、各計算要素ａから複数の方向に直線を伸ばし、各直線が物質表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定する（ステップＳ２４）。計算要素ａから直線を伸ばす方向は、計算要素ａにおける数列θlocal(ｍ),φlocal(ｏ)により決定する。すなわち、ステップＳ２４では、計算要素ａからθlocal(ｍ),φlocal(ｏ)の方向に直線を伸ばす。よって、各計算要素ａからは、Ｍ×Ｏ本の直線を伸ばすこととなる。ステップＳ２４の処理は、Ｎ個の計算要素ａの各々について行われる。ステップＳ２４の処理を行うブロックは、本開示の判定部の例である。
【００６４】
なお、ステップＳ２４では、鏡面境界条件や周期境界条件を考慮して可視判定を行ってもよい。図１３、図１４はそれぞれ、鏡面境界条件、周期境界条件について説明するための模式図である。このような判定を行っておくと、境界条件を取り込んだフラックス計算をローコストで行うことが可能となる。
【００６５】
以上のように、ステップＳ２４では、複数の計算要素ａからの各直線が、物質表面にぶつかるか否かと、どの計算要素にぶつかるかとを判定する。物質表面にぶつかった直線に関しては、ステップＳ２５の処理を行い、物質表面にぶつからなかった直線に関しては、ステップＳ２６の処理を行う。
【００６６】
ステップＳ２５では、ある計算要素ａからのいずれかの直線が計算要素Ｂにぶつかった場合には、その計算要素ａを、計算要素Ｂの可視計算要素としてカウントする。一方、ある計算要素ａからのいずれの直線も計算要素Ｂにぶつからなかった場合には、その計算要素ａは、計算要素Ｂの可視計算要素としてカウントしない。このような処理をすべての計算要素ａについて行うことで、計算要素Ｂから見ることのできるすべての計算要素ａを特定することができる。なお、この処理は、計算要素Ｂのみに限らず、Ｎ個のすべての計算要素について同様に行う。
【００６７】
一方、ステップＳ２６では、ある計算要素ａからのある直線が物質表面にぶつからなかった場合（すなわち、気層に到達した場合）には、その直線の方向を、計算要素ａの気層可視方向としてカウントする。このような処理をすべての直線について行うことで、気層から各計算要素ａに直接的に反応種が到達する方向をすべて特定することができる。この特定結果は、直接フラックスの計算に利用することができる。例えば、計算要素Ｂにおける直接フラックスの計算には、計算要素Ｂについての気層可視方向のカウント結果が利用される。
【００６８】
図８のフローでは次に、ステップＳ２６のカウント結果を利用して、計算要素Ｂにおける直接フラックスΓB,directを計算する（ステップＳ２８）。直接フラックスΓB,directは、数列θlocal(ｍ),φlocal(ｏ)により式（１５）を離散化することで、以下の式（２３）のように表される。
【数２３】
2014135339.tif 000025
【００６９】
ただし、θBlocal(ｍ),φBlocal(ｏ)は、計算要素Ｂにおける数列θlocal(ｍ),φlocal(ｏ)を表す。式（２３）のη(θBlocal,φBlocal)は、計算要素Ｂの気層可視方向ではη＝１となり、その他の方向ではη＝０となる。よって、式（２３）は、ステップＳ２６でカウントされた計算要素Ｂの気層可視方向を利用することで計算することができる。
【００７０】
図８のフローでは次に、ステップＳ２５のカウント結果を利用して、計算要素ａ、Ｂ間の可視係数ν(ａ,Ｂ)と形態係数ｇ(ａ,Ｂ)を計算する（ステップＳ２９）。形態係数ｇ(ａ,Ｂ)は、数列θBlocal(ｍ),φBlocal(ｏ)により以下の式（２４）のように表すことができる。
【数２４】
2014135339.tif 000026
【００７１】
ただし、κ(θBlocal,φBlocal,ａ)は、計算要素ＢからθBlocal,φBlocalの方向に計算要素ａが見えるか否かの可視判定の結果を示しており、計算要素可視判定値と呼ぶことにする。計算要素ＢからθBlocal,φBlocalの方向に計算要素ａが見える場合にはκ(θBlocal,φBlocal,ａ)＝１となり、見えない場合にはκ(θBlocal,φBlocal,ａ)＝０となる。よって、式（２４）は、ステップＳ２５にて計算要素ａが計算要素Ｂの可視計算要素としてカウントされたか否かを参酌することで計算することができる。
【００７２】
計算要素可視判定値κを２次元、３次元で計算する例を、それぞれ図１５、図１６に示す。図１５、図１６はそれぞれ、２次元、３次元における計算要素可視判定値κについて説明するための模式図である。
【００７３】
なお、可視係数ν(ａ,Ｂ)は、式（２４）によるｇ(ａ,Ｂ)の計算結果から算出することができる。具体的には、ｇ(ａ,Ｂ)＝０の場合にはν(ａ,Ｂ)＝０となり、ｇ(ａ,Ｂ)＞０の場合にはν(ａ,Ｂ)＝１となる。
【００７４】
以上のように、ステップＳ２８、Ｓ２９では、ステップＳ２４における判定の結果に基づいて、直接フラックスΓB,direct、可視係数ν(ａ,Ｂ)、および形態係数ｇ(ａ,Ｂ)を計算する。ステップＳ２８、Ｓ２９の処理を行うブロックは、本開示の計算部の例である。なお、ステップＳ２８では、中性種の直接フラックスΓB,ne-directと、イオン種の直接フラックスΓB,ion-directの両方が算出される。
【００７５】
図８のフローによるΓB,ne-direct、ΓB,ion-direct、ν(ａ,Ｂ)、ｇ(ａ,Ｂ)の計算結果は、図１や図４のフローにおいて、総フラックスΓB,ion、ΓB,ne、表面成長速度Ｆi、レベルセット関数ψtなどの計算に用いられる。これらを計算するブロックも、本開示の計算部の例である。
【００７６】
（３）ステップＳ１４の詳細
次に、図１７を参照し、ステップＳ１４の詳細について説明する。
【００７７】
図１７は、図４のステップＳ１４の詳細を示したフローチャート図である。また、図１８〜図２０は、図１７の処理を説明するための模式図である。
【００７８】
ステップＳ１４では、イオン種の反射を取り扱うための反射形態係数ｇionRと、イオン種によるスパッタに起因する中性種の発生を取り扱うためのスパッタ形態係数ｇionSとを計算する。これらの形態係数ｇionR、ｇionSの計算方法は、図８における形態係数ｇの計算方法と概ね同じであるが、以下の２点において形態係数ｇの計算方法と異なっている。
【００７９】
第１に、形態係数ｇionR、ｇionSを計算する際には、図８のステップＳ２４において、直線を伸ばす方向のカットオフを行う（図１９参照）。
【００８０】
図６、図７を参照して説明したように、イオン種は直進性が強く、全方向に反射される訳ではない。そこで、反射形態係数ｇionRの計算用にステップＳ２４を行う際には、イオン種の反射方向のカットオフ角度θcut（図６、図７参照）を設定し、直線を伸ばす方向をカットオフ角度θcutの範囲内に制限する。
【００８１】
また、イオン種の反射と同様に、イオン種によるスパッタに起因する中性種の発生も異方性が強いことが知られている。そこで、スパッタ形態係数ｇionSの計算用にステップＳ２４を行う際には、反射形態係数ｇionRの場合と同様に、中性種の発生方向のカットオフ角度を設定し、直線を伸ばす方向をカットオフ角度の範囲内に制限する。
【００８２】
なお、スパッタ形態係数ｇionS用のカットオフ角度は、反射形態係数ｇionR用のカットオフ角度と同じ値に設定してもよいし、反射形態係数ｇionR用のカットオフ角度とは異なる値に設定してもよい。
【００８３】
図１９では、計算要素ａから直線を伸ばす方向が、カットオフ角度の範囲内に制限されている。その結果、直線がぶつかった計算要素Ｂ1〜Ｂ5は、いずれもカットオフ角度の範囲内に位置している。なお、直線を伸ばす方向は、カットオフ角度以下に制限してもよいし、カットオフ角度未満に制限してもよい。
【００８４】
カットオフ角度θcutは、様々な定義をすることが可能である。例えば、イオン種の入射角度分布を、以下の式（２５）のように定義する場合を想定する。
【数２５】
2014135339.tif 000027
【００８５】
式（２５）は、式（１５）に式（１７）を代入したものに等しいことに留意されたい。
式（２５）を用いる場合、カットオフ角度θcutは、式（２５）における積分対象の式の値が小さくなる方向θlocal、φlocalをカットオフするように設定することが望ましい。この場合、この積分対象の式が計算要素数（計算要素の個数）Ｎに依存しているため、カットオフ角度θcutも、以下の式（２６）のように、計算要素数Ｎに依存するように設定される。
【数２６】
2014135339.tif 000028
【００８６】
式（２６）において、カットオフ角度θcutは、計算要素数Ｎの関数であり、計算要素数Ｎに依存している。この場合には、例えば、Ｎ＝１０のときにはθcut＝３０度、Ｎ＝１００のときにはθcut＝１０度、Ｎ＝１０００のときにはθcut＝３度のように、カットオフ角度θcutが計算要素数Ｎに応じて変化することとなる。
【００８７】
第２に、形態係数ｇionR、ｇionSを計算する際に、計算要素ａからの直線が計算要素Ｂにぶつかった場合には、計算要素Ｂの周囲の計算要素Ｃ、Ｃ’についてもステップＳ２４の処理を行う（図１８、図２０参照）。計算要素Ｃは、計算要素Ｂに直接的に隣接する計算要素であり、計算要素Ｃ’は、計算要素Ｂに計算要素Ｃを介して間接的に隣接する計算要素である。計算要素ａ、計算要素Ｂ、計算要素Ｃ、Ｃ’はそれぞれ、本開示の第１、第２、第３の計算要素の例である。
【００８８】
具体的にはまず、計算要素ａから、計算要素Ｂに直接的に隣接する計算要素Ｃに向けて新たな直線を伸ばし、この直線が他の計算要素を介さずに計算要素Ｃにぶつかるか否かを判定する。また、計算要素Ｃが、計算要素ａのカットオフ角度θcutの範囲内に位置するか否かも判定する。これにより、計算要素Ｃについても、計算要素Ｂと同様にステップＳ２４の判定処理が行われたこととなる。
【００８９】
また、計算要素ａからの直線が計算要素Ｃにぶつかり、かつ、計算要素Ｃが計算要素ａのカットオフ角度θcutの範囲内に位置する場合（すなわち、計算要素Ｃについて肯定の判定結果が得られた場合）には、計算要素Ｃに直接的に隣接する計算要素Ｃ’についてもステップＳ２４の処理を行う。
【００９０】
本実施形態では、このような処理を、判定対象の計算要素がなくなるまで繰り返す。すなわち、本実施形態では、第３の計算要素として、第２の計算要素Ｂに直接的に隣接する計算要素と、第２の計算要素Ｂに、肯定の判定結果が得られた計算要素を介して間接的に隣接する計算要素とを選択し、選択可能な第３の計算要素がなくなるまで判定処理を繰り返す。
【００９１】
図２０では、図１９の計算要素Ｂ1〜Ｂ5に加え、計算要素Ｃ1〜Ｃ12についてもステップＳ２４の処理が行われている。計算要素ａからの直線は、計算要素Ｃ1〜Ｃ5、Ｃ8、Ｃ10、Ｃ11にぶつかるが、そのうち、計算要素Ｃ1はカットオフ角度の範囲外に位置しており、形態係数ｇionR、ｇionSの計算から除外される。また、計算要素Ｃ6、Ｃ7、Ｃ9、Ｃ12は、計算要素ａから見て他の計算要素の背後にあるため、計算要素ａからの直線がぶつからず、形態係数ｇionR、ｇionSの計算から除外される。
【００９２】
なお、形態係数ｇionR、ｇionSの計算時には図２０の処理が行われるため、形態係数ｇの計算時よりステップＳ２１の分割数（直線を伸ばす本数）Ｍ×Ｏを減らしても、十分な計算精度を得ることが可能である。これにより、本実施形態では、形態係数ｇionR、ｇionSの計算時におけるステップＳ２４の計算時間を、形態係数ｇの計算時よりも短縮することが可能となる。また、本実施形態では、形態係数ｇionR、ｇionSの計算時に、図２０の処理により計算要素Ｂだけでなく計算要素Ｂの周囲の計算要素でも判定処理を行うことで、判定処理をより細かく行い、計算誤差を低減することが可能となる。
【００９３】
図１７および図１８は、図２０の処理の詳細を示している。以下、図１７および図１８を参照し、図２０の処理の詳細について説明する。図１７および図１８は、図２０の処理をシード・フィルアルゴリズムを用いて実行する一例を示している。
【００９４】
図１７のステップＳ３１〜Ｓ３３では、計算要素ａからカットオフ角度の範囲内の複数の方向に直線を伸ばし、各直線がどの計算要素にぶつかるかを判定している。すなわち、ステップＳ３１〜Ｓ３３は、図８のステップＳ２１〜Ｓ２４に相当している。
【００９５】
ステップＳ３３では、計算要素ａからの直線が計算要素Ｂにぶつかるか否かを判定している。この処理の様子を、図１８(ａ)に示す。図１８(ａ)の四角形は、計算要素Ｂとその周囲の計算要素を示す。各四角形内の数値は、各計算要素に立てられたフラグを示す。
【００９６】
フラグ「０」は、初期値に相当する。また、フラグ「１」が立てられた計算要素は、計算要素ａからの直線がその計算要素にぶつかり、かつ、その計算要素がカットオフ角度の範囲内に位置することを意味する。また、フラグ「２」が立てられた計算要素は、計算要素ａからの直線がその計算要素にぶつからないか、その計算要素がカットオフ角度の範囲外に位置することを意味する。
【００９７】
計算要素ａからの直線が計算要素Ｂにぶつかった場合には、計算要素Ｂにフラグ「１」が立てられる（ステップＳ３５）。一方、計算要素ａからの直線が計算要素Ｂにぶつからなかった場合には、計算要素Ｂにフラグ「２」が立てられる（ステップＳ４０）。図８(ａ)は、計算要素Ｂにフラグ「１」が立てられた様子を示している。
【００９８】
なお、ステップＳ３４では、計算要素Ｂが計算要素ａのカットオフ角度の範囲内に位置するか否かを判定している。しかしながら、ステップＳ３１〜Ｓ３２ではカットオフ角度の範囲内の方向のみに直線を伸ばすため、計算要素Ｂは原則的にカットオフ角度の範囲内に位置している。このステップＳ３４は、後に計算要素Ｂの周囲の計算要素についてカットオフ判定を行う際に意味が出てくる。なお、図１７に示すＸmax、Ｙmax、Ｘmin、Ｙminはそれぞれ、計算要素Ｂの最大Ｘ座標、最大Ｙ座標、最小Ｘ座標、最小Ｙ座標を示す。
【００９９】
計算要素ａからの直線が計算要素Ｂにぶつかった場合には、図１８(ａ)に示すように、計算要素Ｂに直接的に隣接する各計算要素Ｃについても、計算要素Ｂと同様の処理を行う（ステップＳ３６〜Ｓ３９）。
【０１００】
具体的には、計算要素ａから計算要素Ｃに向けて新たな直線を伸ばし、この直線が計算要素Ｃにぶつかるか否かを判定する（ステップＳ３３）。また、計算要素Ｃが、計算要素ａのカットオフ角度の範囲内に位置するか否かも判定する（ステップＳ３４）。
【０１０１】
そして、計算要素ａからの直線がある計算要素Ｃにぶつかり、かつ、その計算要素Ｃがカットオフ角度の範囲内に位置する場合には、その計算要素Ｃにフラグ「１」が立てられる（ステップＳ３５）。一方、計算要素ａからの直線がその計算要素Ｃにぶつからなかった場合や、その計算要素Ｃがカットオフ角度の範囲外に位置する場合には、その計算要素Ｃにフラグ「２」が立てられる（ステップＳ４０）。図８(ｂ)は、各計算要素Ｃにフラグ「１」または「２」が立てられた様子を示している。
【０１０２】
ある計算要素Ｃにフラグ「１」が立てられた場合には、図１８(ｃ)に示すように、その計算要素Ｃに直接的に隣接する各計算要素Ｃ’についても、計算要素Ｂと同様の処理を行う（ステップＳ３６〜Ｓ３９）。ただし、すでにフラグ「１」または「２」が立っている計算要素Ｃ’については、この処理は不要である。
【０１０３】
本実施形態では、ステップＳ３６〜Ｓ３９の処理を、判定対象の計算要素がなくなるまで繰り返す。具体的には、図１８(ｄ)に示すように、フラグ「１」が立っている計算要素が、フラグ「２」が立っている計算要素に囲まれた状態になるまで、ステップＳ３６〜Ｓ３９の処理が繰り返される。
【０１０４】
図１８(ｄ)の符号Ｒは、フラグ「１」が立っている計算要素からなる領域を示す。本実施形態では、この領域Ｒ内に含まれる計算要素を、計算要素ａからの直線がぶつかり、かつ、計算要素ａのカットオフ角度の範囲内に位置する計算要素としてカウントする。
【０１０５】
なお、本実施形態では、ステップＳ１２〜Ｓ１４において、各計算要素に固有のローカル座標系を使用したが、代わりに、すべての計算要素に共通のグローバル座標系を使用してもよい。
【０１０６】
（４）第１実施形態における計算時間と計算誤差
次に、以上の説明を踏まえて、第１実施形態における計算時間と計算誤差について説明する。
【０１０７】
従来の方法では、任意の計算要素Ｂの直接フラックスΓB,ne-direct、ΓB,ion-directを計算するのに、計算要素数Ｎに比例する時間がかかる。理由は、計算要素Ｂに関するループ計算をＮ回繰り返し行うからである。また、従来の方法では、任意の計算要素ａ、Ｂ間の可視係数ν(ａ,Ｂ)や形態係数ｇ(ａ,Ｂ)、ｇionR(ａ,Ｂ)、ｇionS(ａ,Ｂ)を計算するのに、Ｎ2に比例する時間がかかる。理由は、計算要素ａに関するループ計算と計算要素Ｂに関するループ計算を、それぞれＮ回繰り返し行うからである。可視係数や形態係数の計算時間は、鏡面境界条件や周期境界条件を採用するとさらに長くなる。よって、従来の方法における計算時間の多くは、可視係数や形態係数の計算に費やされる。
【０１０８】
一方、本実施形態では、図８や図１７に示すように、各計算要素ａから複数の方向に直線を伸ばし、各直線が物質表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定し、この判定結果に基づいて直接フラックス、可視係数、形態係数を計算する。そのため、可視係数と形態係数が、直接フラックスと同様に、計算要素ａに関するループ計算をＮ回繰り返すことで計算される（ステップＳ２２、Ｓ３０を参照）。よって、本実施形態によれば、直接フラックス、可視係数、形態係数の計算時間を、計算要素数Ｎに比例する時間に抑えることができる。
【０１０９】
また、本実施形態では、イオン種と中性種の特性の違いを考慮し、反射形態係数ｇionRやスパッタ形態係数ｇionSの計算用にステップＳ２４の判定処理を行う際には、イオン種の反射方向や、イオン種によるスパッタに起因する中性種の発生方向のカットオフ角度を設定し、直線を伸ばす方向をカットオフ角度の範囲内に制限する。さらに、本実施形態では、上記の判定処理を、これらの直線がぶつかった計算要素の周囲の計算要素についても繰り返し適用する。よって、本実施形態によれば、計算の無駄を減らすことで計算時間を短縮することや、無駄な計算の代わりに有用な計算に時間を割くことで計算誤差を低減することが可能となる。なお、中性種の発生量が少ないなど、中性種の発生を無視できる場合には、式（４）におけるスパッタ形態係数ｇionSを含む項を無視して形状シミュレーションを行ってもよい。
【０１１０】
また、本実施形態によるｇ、ｇionR、ｇionSの計算によれば、Ｎ2回のループ計算でｇ、ｇionR、ｇionSを計算する従来の方法に比べて、ｇ行列、ｇionR行列、ｇionS行列（さらにはν行列）における０要素の個数が多くなる傾向にある。本実施形態では、各計算要素から複数の方向に直線を伸ばし、各直線が物質表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定することで、形態係数を計算する。その結果、本実施形態によれば、すべての計算要素の組合せ間でループ計算を行う場合と比べて、形態係数が０となる割合が大幅に増え、ｇ行列、ｇionR行列、ｇionS行列の各々の全行列要素に占める０要素の割合は、半数以上（具体的には、多くの場合８割程度）となる。この場合、式（１３）の行列Ａneの非対角要素は半数以上が０となり、式（１０）の行列方程式はシンプルな形となる（同様に、イオン種に関する行列Ａionの非対角要素も半数以上が０となり、行列Ａionを含む行列方程式もシンプルな形となる）。この結果として、本実施形態によれば、計算時間とメモリ使用量を大幅に減らすことが可能となる。
【０１１１】
よって、本実施形態によれば、これらの行列方程式を繰り返し解きながら化学反応計算を行う場合に、これらの０要素に着目した計算アルゴリズムを採用することで、計算時間をさらに短縮することが可能となる。さらには、ＣＲＳなどの疎行列保持アルゴリズムを採用することで、０要素が多くなるほど省メモリ化を図ることも可能となる。この場合には、図４のステップＳ１７において、Ｓi(Γi)、Ｒi(Γi)、Ｐi(Γi)が収束するまで行列方程式を繰り返し解くことになるが、０要素が多いことで行列方程式を１回解くのに要する計算時間が短縮されるため、ステップＳ１７のトータルの計算時間が大幅に短縮される。
【０１１２】
（５）第１実施形態の効果
最後に、第１実施形態の効果について説明する。
【０１１３】
以上のように、本実施形態では、各計算要素から複数の方向に直線を伸ばし、各直線が物質表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定し、この判定結果に基づいて直接フラックスと形態係数を計算する。さらには、この判定結果に基づいて可視係数を計算する。
【０１１４】
よって、本実施形態によれば、直接フラックスや形態係数の計算時間を、計算要素数に比例する時間に抑えることができる。よって、本実施形態によれば、間接フラックスの計算時間に影響を与える形態係数の計算時間を短縮することで、物質表面に直接的、間接的に到達する反応種を考慮した形状シミュレーションを高速化することが可能となる。
【０１１５】
この効果の具体例については、図２１〜図２４を参照して後述する。
【０１１６】
また、本実施形態では、イオン種と中性種の特性の違いを考慮し、反射形態係数やスパッタ形態係数の計算用に上記の判定処理を行う場合には、イオン種の反射方向や、イオン種によるスパッタに起因する中性種の発生方向のカットオフ角度を設定し、直線を伸ばす方向をカットオフ角度の範囲内に制限する。さらには、上記の判定処理を、これらの直線がぶつかった計算要素の周囲の計算要素についても繰り返し適用する。
【０１１７】
よって、本実施形態によれば、計算の無駄を減らして計算時間を短縮することや、無駄な計算の代わりに有用な計算に時間を割いて計算誤差を低減することなどが可能となり、形状シミュレーションを高速化および高精度化することが可能となる。
【０１１８】
この効果の具体例については、図２５〜図２８を参照して後述する。
【０１１９】
（５．１）図２１〜図２４の説明
図２１、図２２はそれぞれ、比較例、第１実施形態における計算時間の例を示したグラフである。ただし、図２１、図２２では、カットオフ角度の設定前における第１実施形態の効果を検証するために、反射形態係数ｇionRとスパッタ形態係数ｇionSを無視して計算を行った（図２３、図２４も同様）。そして、比較例では、従来の方法を用いて直接フラックス、可視係数、形態係数（ｇ）を計算した。図２１、図２２は、図２に示す構造を初期構造とした場合の直接フラックスの計算時間、可視計算（可視係数と形態係数の計算）の計算時間、化学反応収束計算の計算時間、全計算時間の合計を示している。
【０１２０】
図２１、図２２に示すように、第１実施形態によれば、全計算時間を比較例よりも短縮することが可能となる。これらの比較結果を図２３に示す。図２３は、計算要素数が４万の場合の第１実施形態と比較例の計算時間を比較したグラフである。
【０１２１】
図２４は、第１実施形態と比較例におけるθ分割数と計算誤差との関係を示したグラフである。ただし、図２４における計算では、ローカル座標系を使用した。また、ｎ＝１とｎ＝１０００のグラフは、比較例による計算結果を示す。図２４に示すように、第１実施形態によれば、比較例よりも計算誤差を抑制できることが分かる。
【０１２２】
（５．２）図２５〜図２８の説明
図２５は、イオン種のみを取り扱った場合の第１実施形態と比較例の計算時間を比較したグラフである。また、図２６は、イオン種と中性種を取り扱った場合の第１実施形態と比較例の計算時間を比較したグラフである。図２５では、１種類のイオンを取り扱い、図２６では、１種類のイオンと１種類の中性粒子とを取り扱った。また、図２５、図２６では、第１実施形態において、カットオフ角度の設定や図１８の処理を行った（図２７、図２８も同様）。図２５、図２６は、図２に示す構造を初期構造とした場合の１ステップ当たりの全計算時間を示している。
【０１２３】
図２５、図２６に示すように、第１実施形態によれば、図２１、図２２の場合よりもさらに顕著に、全計算時間を比較例よりも短縮することが可能となる。これらの比較結果の詳細を図２７、図２８に示す。
【０１２４】
図２７は、イオン種のみを取り扱った場合の第１実施形態と比較例の計算時間を項目別にに比較したグラフであり、図２５に対応している。また、図２８は、イオン種と中性種を取り扱った場合の第１実施形態と比較例の計算時間を項目別に比較したグラフであり、図２６に対応している。図２７、図２８における計算要素数は、約４万である。図２７、図２８に示す例では、第１実施形態における間接フラックスの計算時間を、比較例に比べて大幅に短縮することができた。
【０１２５】
なお、第１実施形態の形状シミュレーション方法は、どのような情報処理装置を用いて実行してもよい。第２実施形態では、このような情報処理装置の一例である形状シミュレーション装置について説明する。
【０１２６】
（第２実施形態）
図２９は、第２実施形態の形状シミュレーション装置の構成を示す外観図である。
【０１２７】
図２９の形状シミュレーション装置は、制御部１１と、表示部１２と、入力部１３とを備えている。
【０１２８】
制御部１１は、形状シミュレーション装置の動作を制御するモジュールである。制御部１１は、例えば、第１実施形態の形状シミュレーション方法を実行する。制御部１１の詳細については後述する。
【０１２９】
表示部１２は、液晶モニタなどの表示デバイスを有している。表示部１２は、例えば、形状シミュレーション用の設定情報の入力画面や、形状シミュレーションの計算結果などを表示する。
【０１３０】
入力部１３は、キーボード１３ａやマウス１３ｂなどの入力デバイスを有している。入力部１３は、例えば、形状シミュレーション用の設定情報の入力用に使用される。設定情報の例としては、計算式に関する情報、実験値や予測値に関する情報、物質の構造に関する情報、フラックスに関する情報、形状シミュレーションの条件や手順に関する指示情報などが挙げられる。
【０１３１】
図３０は、図２９の制御部１１の構成を示すブロック図である。
【０１３２】
制御部１１は、ＣＰＵ（Central Processing Unit）２１と、ＲＯＭ（Read Only Memory）２２と、ＲＡＭ（Random Access Memory）２３と、ＨＤＤ（Hard Disk Drive）２４と、ＣＤ（Compact Disc）ドライブなどのメモリドライブ２５と、メモリポートやメモリスロットなどのメモリＩ／Ｆ（interface）２６とを備えている。
【０１３３】
本実施形態では、第１実施形態の形状シミュレーション方法用のプログラムである形状シミュレーションプログラムが、ＲＯＭ２２またはＨＤＤ２４内に格納されている。入力部１３から所定の指示情報が入力されると、ＣＰＵ２１は、ＲＯＭ２２またはＨＤＤ２４からプログラムを読み出し、読み出したプログラムをＲＡＭ２３に展開し、このプログラムにより形状シミュレーションを実行する。この処理の際に生じる各種データは、ＲＡＭ２３内に保持される。
【０１３４】
なお、本実施形態では、形状シミュレーションプログラムを記録したコンピュータ読み取り可能な記録媒体を用意し、この記録媒体からＲＯＭ２２やＨＤＤ２４内に形状シミュレーションプログラムをインストールしてもよい。このような記録媒体の例としては、ＣＤ−ＲＯＭやＤＶＤ（Digital Versatile Disk）−ＲＯＭなどが挙げられる。
【０１３５】
また、本実施形態では、形状シミュレーションプログラムを、インターネットなどのネットワーク経由でダウンロードすることで、ＲＯＭ２２やＨＤＤ２４内にインストールしてもよい。
【０１３６】
以上のように、本実施形態によれば、第１実施形態の形状シミュレーション方法を実行するための形状シミュレーション装置や形状シミュレーションプログラムを提供することが可能となる。
【０１３７】
なお、第１、第２実施形態では、形状シミュレーションの適用対象の例として半導体デバイスを取り上げたが、この形状シミュレーションは半導体デバイス以外のデバイスにも適用可能である。このようなデバイスの例としては、ＭＥＭＳ（Micro Electro Mechanical Systems）デバイスやディスプレイデバイスなどが挙げられる。
【０１３８】
以上、いくつかの実施形態を説明したが、これらの実施形態は、例としてのみ提示したものであり、発明の範囲を限定することを意図したものではない。本明細書で説明した新規な装置、方法、およびプログラムは、その他の様々な形態で実施することができる。また、本明細書で説明した装置、方法、およびプログラムの形態に対し、発明の要旨を逸脱しない範囲内で、種々の省略、置換、変更を行うことができる。添付の特許請求の範囲およびこれに均等な範囲は、発明の範囲や要旨に含まれるこのような形態や変形例を含むように意図されている。
【符号の説明】
【０１３９】
１：シリコン基板、２：シリコン窒化膜、３：シリコン酸化膜、４：貫通孔、
１１：制御部、１２：表示部、１３：入力部、
２１：ＣＰＵ、２２：ＲＯＭ、２３：ＲＡＭ、２４：ＨＤＤ、
２５：メモリドライブ、２６：メモリＩ／Ｆ
【図１】
2014135339.tif 000029
【図２】
2014135339.tif 000030
【図３】
2014135339.tif 000031
【図４】
2014135339.tif 000032
【図５】
2014135339.tif 000033
【図６】
2014135339.tif 000034
【図７】
2014135339.tif 000035
【図８】
2014135339.tif 000036
【図９】
2014135339.tif 000037
【図１０】
2014135339.tif 000038
【図１１】
2014135339.tif 000039
【図１２】
2014135339.tif 000040
【図１３】
2014135339.tif 000041
【図１４】
2014135339.tif 000042
【図１５】
2014135339.tif 000043
【図１６】
2014135339.tif 000044
【図１７】
2014135339.tif 000045
【図１８】
2014135339.tif 000046
【図１９】
2014135339.tif 000047
【図２０】
2014135339.tif 000048
【図２１】
2014135339.tif 000049
【図２２】
2014135339.tif 000050
【図２３】
2014135339.tif 000051
【図２４】
2014135339.tif 000052
【図２５】
2014135339.tif 000053
【図２６】
2014135339.tif 000054
【図２７】
2014135339.tif 000055
【図２８】
2014135339.tif 000056
【図２９】
2014135339.tif 000057
【図３０】
2014135339.tif 000058