JP6034700B2 - 形状シミュレーション装置、形状シミュレーション方法、および形状シミュレーションプログラム - Google Patents

Description

本発明の実施形態は、形状シミュレーション装置、形状シミュレーション方法、および形状シミュレーションプログラムに関する。

ＣＶＤ（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は、以下の式（１）から計算できる。

ただし、∇はベクトル微分演算子を表し、｜∇ψ_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}の合計との和で表される。

同様に、計算要素Ｂにおけるイオン種の総フラックスΓ_B,ionは、以下の式（３）のように、計算要素Ｂにおけるイオン種の直接フラックスΓ_B,ion-directと、任意の計算要素ａからのイオン種の間接フラックスΓ_{aB,ion-indirect}の合計との和で表される。

ここで、間接フラックスΓ_{aB,ne-indirect}、Γ_{aB,ion-indirect}は、例えば以下の式（４）、式（５）で表すことができる。

Ｓ_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はそれぞれ、以下の式（６）、式（７）で表すことができる。

図４のフローでは次に、任意の計算要素における直接フラックスΓ_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を計算する（ステップＳ１６）。

次に、ステップＳ１５とステップＳ１６の処理を、付着確率Ｓ_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を、それぞれ可視係数行列、形態係数行列と呼ぶ。また、任意の計算要素におけるフラックスは、まとめてＮ行ベクトルで表すことができる。ベクトルの形で表されたフラックスを、フラックスベクトルと呼ぶ。

この場合、式（８）は、以下の式（１０）のように行列方程式で表現することができる。

ただし、Ｉ、Ｊは処理対象の計算要素の個数を表し、例えばＩ＝Ｊ＝Ｎである。また、式（１３）では、表記スペースの都合上、付着確率Ｓ_j(Γ_j,ion,Γ_j,ne)をＳ_j(Γ_j)と略記すると共に、スパッタ確率Ｐ_i(Γ_i,ion,Γ_i,ne)を含む項の記載を省略している。

行列方程式（１０）は、どのような解法で解いてもよい。解法の例としては、反復法（ガウスザイデル法、ＳＯＲ法、ヤコビ法、共役勾配法など）や、直接法（ガウスの消去法、ＬＵ分解法、コレスキー分解法など）が挙げられる。行列方程式（１０）を解く際に、行列Ａ_neが疎行列である場合には、ＣＲＳなどの保存方法を用いた上で疎行列に適したルーチンを使用することで、計算処理の省メモリ化や高速化を図ってもよい。

なお、式（９）も、式（８）と同様に行列方程式で表現可能である。本実施形態では、これら２つの行列方程式を上記の解法で解くことができる。

図４のフローでは次に、総フラックスΓ_i,ion、Γ_i,neから、各計算要素ｉにおける局所的な表面成長速度Ｆ_iを計算する（ステップＳ１８）。例えば、Ｋ種類の反応種を使用する場合においては、表面成長速度Ｆ_iは、Ｋ個の局所的な総フラックスΓ_1,i〜Γ_K,iに依存した以下の式（１４）のような形でモデリングされる。

ただし、ｋは１≦ｋ≦Ｋを満たす任意の実数である。なお、Ｋ種類の反応種は、中性種とイオン種のいずれか一方のみを含んでいてもよいし、中性種とイオン種の両方を含んでいてもよい。以上のようにして、ステップＳ３の処理が終了する。

（２）ステップＳ１２、Ｓ１３の詳細
次に、図８を参照し、ステップＳ１２、Ｓ１３の詳細について説明する。

ステップＳ１２、Ｓ１３では、中性種の直接フラックスΓ_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は、以下の式（１５）により計算することができる。

ただし、η(θ_local,φ_local)は、計算要素Ｂからθ_local,φ_localの方向に直線を伸ばした場合の可視判定の結果を示しており、可視判定値と呼ぶことにする。図１０は、可視判定値ηについて説明するための模式図である。図１０に示すように、上記直線が物質表面とぶつかる場合には、η＝０となり、ぶつからない場合にはη＝１となる。なお、図１０に示すように、物質表面の片側の方向にだけ直線を伸ばす場合には、式（１５）におけるθ_localの積分範囲は、０からπではなく、０からπ／２にしてもよい。

可視判定値ηと可視係数νとの違いについては、図１１を参照されたい。図１１は、可視係数νについて説明するための模式図である。ν(ａ,Ｂ)は、計算要素ａと計算要素Ｂが互いに見えるか否かを示す。計算要素ａ、Ｂを結ぶ直線が、計算要素ａ、Ｂ間において物質表面とぶつかる場合には、ν＝０となり、ぶつからない場合にはν＝１となる。前者の例としては計算要素ｄを、後者の例としては計算要素ｃを参照されたい。

また、ｆ_flatは、平坦面での直接フラックスを示し、入力値として事前に与えられる。さらに、Normは、以下の式（１６）で与えられる規格化定数を示す。さらに、ｆ(θ_local)は、直接フラックスの面積素片の係数を示し、例えば以下の式（１７）で与えられる。

ただし、θ_inは、図１２に示すような入射角度である。図１２は、入射角度θ_inについて説明するための模式図である。入射角度θ_inは、法線ベクトル方向とθ_local,φ_localの方向との間の角度に相当する。よって、ローカル座標系（ｒ_local,θ_local,φ_local）を用いる場合には、θ_in＝θ_localが成り立つ。

以下、図８のフローについて具体的に説明する。

まず、極角θ_localの数列θ_local(ｍ)の値（ｍ＝０、１、…、Ｍ−１）と、偏角φ_localの数列φ_local(ｏ)の値（ｏ＝０、１、…、Ｏ−１）とを計算する（ステップＳ２１）。これは、０からπまでの極角θ_localの範囲をＭ個の領域に分割し、０から２πまでの偏角φ_localの範囲をＯ個の領域に分割することに相当する。後述するように、式（１５）の積分計算は、これらの数列θ_local(ｍ),φ_local(ｏ)を用いて離散化される。

式（１５）の直接フラックスΓ_B,directの計算に、式（１７）に示す面積素片係数を使用する場合には、例えば、以下の式（１８）、式（１９）のような数列θ_local(ｍ),φ_local(ｏ)を用意する。

ここで、数列∂(ｍ)は、以下の式（２０）で与えられる。

式（１８）のθ_local(ｍ)は、ｆ(θ_local)｜sinθ_local｜をθ_local＝０からθ_local＝θ_local(ｍ)まで積分した場合に、積分結果が∂(ｍ)となる角度を表す。この定義から式（２１）の関係が成り立ち、式（２１）から式（２２）が導出され、式（２２）を変形することで式（１８）が得られる。

以上のように、ステップＳ２１では、０からπまでの極角θ_localの範囲を非等間隔に分割し、０から２πまでの偏角φ_localの範囲を等間隔に分割する。なお、本実施形態では、極角θ_localの範囲だけでなく、偏角φ_localの範囲も非等間隔に分割してもよい。また、極角θ_localの積分範囲を０からπ／２とする場合には、０からπまでではなく、０からπ／２までの極角θ_localの範囲をＭ個の領域に分割するようにしてもよい。

次に、各計算要素ａから複数の方向に直線を伸ばし、各直線が物質表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定する（ステップＳ２４）。計算要素ａから直線を伸ばす方向は、計算要素ａにおける数列θ_local(ｍ),φ_local(ｏ)により決定する。すなわち、ステップＳ２４では、計算要素ａからθ_local(ｍ),φ_local(ｏ)の方向に直線を伸ばす。よって、各計算要素ａからは、Ｍ×Ｏ本の直線を伸ばすこととなる。ステップＳ２４の処理は、Ｎ個の計算要素ａの各々について行われる。ステップＳ２４の処理を行うブロックは、本開示の判定部の例である。

なお、ステップＳ２４では、鏡面境界条件や周期境界条件を考慮して可視判定を行ってもよい。図１３、図１４はそれぞれ、鏡面境界条件、周期境界条件について説明するための模式図である。このような判定を行っておくと、境界条件を取り込んだフラックス計算をローコストで行うことが可能となる。

以上のように、ステップＳ２４では、複数の計算要素ａからの各直線が、物質表面にぶつかるか否かと、どの計算要素にぶつかるかとを判定する。物質表面にぶつかった直線に関しては、ステップＳ２５の処理を行い、物質表面にぶつからなかった直線に関しては、ステップＳ２６の処理を行う。

ステップＳ２５では、ある計算要素ａからのいずれかの直線が計算要素Ｂにぶつかった場合には、その計算要素ａを、計算要素Ｂの可視計算要素としてカウントする。一方、ある計算要素ａからのいずれの直線も計算要素Ｂにぶつからなかった場合には、その計算要素ａは、計算要素Ｂの可視計算要素としてカウントしない。このような処理をすべての計算要素ａについて行うことで、計算要素Ｂから見ることのできるすべての計算要素ａを特定することができる。なお、この処理は、計算要素Ｂのみに限らず、Ｎ個のすべての計算要素について同様に行う。

一方、ステップＳ２６では、ある計算要素ａからのある直線が物質表面にぶつからなかった場合（すなわち、気層に到達した場合）には、その直線の方向を、計算要素ａの気層可視方向としてカウントする。このような処理をすべての直線について行うことで、気層から各計算要素ａに直接的に反応種が到達する方向をすべて特定することができる。この特定結果は、直接フラックスの計算に利用することができる。例えば、計算要素Ｂにおける直接フラックスの計算には、計算要素Ｂについての気層可視方向のカウント結果が利用される。

図８のフローでは次に、ステップＳ２６のカウント結果を利用して、計算要素Ｂにおける直接フラックスΓ_B,directを計算する（ステップＳ２８）。直接フラックスΓ_B,directは、数列θ_local(ｍ),φ_local(ｏ)により式（１５）を離散化することで、以下の式（２３）のように表される。

ただし、θ_Blocal(ｍ),φ_Blocal(ｏ)は、計算要素Ｂにおける数列θ_local(ｍ),φ_local(ｏ)を表す。式（２３）のη(θ_Blocal,φ_Blocal)は、計算要素Ｂの気層可視方向ではη＝１となり、その他の方向ではη＝０となる。よって、式（２３）は、ステップＳ２６でカウントされた計算要素Ｂの気層可視方向を利用することで計算することができる。

図８のフローでは次に、ステップＳ２５のカウント結果を利用して、計算要素ａ、Ｂ間の可視係数ν(ａ,Ｂ)と形態係数ｇ(ａ,Ｂ)を計算する（ステップＳ２９）。形態係数ｇ(ａ,Ｂ)は、数列θ_Blocal(ｍ),φ_Blocal(ｏ)により以下の式（２４）のように表すことができる。

ただし、κ(θ_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用のカットオフ角度とは異なる値に設定してもよい。

図１９では、計算要素ａから直線を伸ばす方向が、カットオフ角度の範囲内に制限されている。その結果、直線がぶつかった計算要素Ｂ₁〜Ｂ₅は、いずれもカットオフ角度の範囲内に位置している。なお、直線を伸ばす方向は、カットオフ角度以下に制限してもよいし、カットオフ角度未満に制限してもよい。

カットオフ角度θ_cutは、様々な定義をすることが可能である。例えば、イオン種の入射角度分布を、以下の式（２５）のように定義する場合を想定する。

式（２５）は、式（１５）に式（１７）を代入したものに等しいことに留意されたい。
式（２５）を用いる場合、カットオフ角度θ_cutは、式（２５）における積分対象の式の値が小さくなる方向θ_local、φ_localをカットオフするように設定することが望ましい。この場合、この積分対象の式が計算要素数（計算要素の個数）Ｎに依存しているため、カットオフ角度θ_cutも、以下の式（２６）のように、計算要素数Ｎに依存するように設定される。

式（２６）において、カットオフ角度θ_cutは、計算要素数Ｎの関数であり、計算要素数Ｎに依存している。この場合には、例えば、Ｎ＝１０のときにはθ_cut＝３０度、Ｎ＝１００のときにはθ_cut＝１０度、Ｎ＝１０００のときにはθ_cut＝３度のように、カットオフ角度θ_cutが計算要素数Ｎに応じて変化することとなる。

第２に、形態係数ｇ_ionR、ｇ_ionSを計算する際に、計算要素ａからの直線が計算要素Ｂにぶつかった場合には、計算要素Ｂの周囲の計算要素Ｃ、Ｃ’についてもステップＳ２４の処理を行う（図１８、図２０参照）。計算要素Ｃは、計算要素Ｂに直接的に隣接する計算要素であり、計算要素Ｃ’は、計算要素Ｂに計算要素Ｃを介して間接的に隣接する計算要素である。計算要素ａ、計算要素Ｂ、計算要素Ｃ、Ｃ’はそれぞれ、本開示の第１、第２、第３の計算要素の例である。

具体的にはまず、計算要素ａから、計算要素Ｂに直接的に隣接する計算要素Ｃに向けて新たな直線を伸ばし、この直線が他の計算要素を介さずに計算要素Ｃにぶつかるか否かを判定する。また、計算要素Ｃが、計算要素ａのカットオフ角度θ_cutの範囲内に位置するか否かも判定する。これにより、計算要素Ｃについても、計算要素Ｂと同様にステップＳ２４の判定処理が行われたこととなる。

また、計算要素ａからの直線が計算要素Ｃにぶつかり、かつ、計算要素Ｃが計算要素ａのカットオフ角度θ_cutの範囲内に位置する場合（すなわち、計算要素Ｃについて肯定の判定結果が得られた場合）には、計算要素Ｃに直接的に隣接する計算要素Ｃ’についてもステップＳ２４の処理を行う。

本実施形態では、このような処理を、判定対象の計算要素がなくなるまで繰り返す。すなわち、本実施形態では、第３の計算要素として、第２の計算要素Ｂに直接的に隣接する計算要素と、第２の計算要素Ｂに、肯定の判定結果が得られた計算要素を介して間接的に隣接する計算要素とを選択し、選択可能な第３の計算要素がなくなるまで判定処理を繰り返す。

図２０では、図１９の計算要素Ｂ₁〜Ｂ₅に加え、計算要素Ｃ₁〜Ｃ₁₂についてもステップＳ２４の処理が行われている。計算要素ａからの直線は、計算要素Ｃ₁〜Ｃ₅、Ｃ₈、Ｃ₁₀、Ｃ₁₁にぶつかるが、そのうち、計算要素Ｃ₁はカットオフ角度の範囲外に位置しており、形態係数ｇ_ionR、ｇ_ionSの計算から除外される。また、計算要素Ｃ₆、Ｃ₇、Ｃ₉、Ｃ₁₂は、計算要素ａから見て他の計算要素の背後にあるため、計算要素ａからの直線がぶつからず、形態係数ｇ_ionR、ｇ_ionSの計算から除外される。

なお、形態係数ｇ_ionR、ｇ_ionSの計算時には図２０の処理が行われるため、形態係数ｇの計算時よりステップＳ２１の分割数（直線を伸ばす本数）Ｍ×Ｏを減らしても、十分な計算精度を得ることが可能である。これにより、本実施形態では、形態係数ｇ_ionR、ｇ_ionSの計算時におけるステップＳ２４の計算時間を、形態係数ｇの計算時よりも短縮することが可能となる。また、本実施形態では、形態係数ｇ_ionR、ｇ_ionSの計算時に、図２０の処理により計算要素Ｂだけでなく計算要素Ｂの周囲の計算要素でも判定処理を行うことで、判定処理をより細かく行い、計算誤差を低減することが可能となる。

図１７および図１８は、図２０の処理の詳細を示している。以下、図１７および図１８を参照し、図２０の処理の詳細について説明する。図１７および図１８は、図２０の処理をシード・フィルアルゴリズムを用いて実行する一例を示している。

図１７のステップＳ３１〜Ｓ３３では、計算要素ａからカットオフ角度の範囲内の複数の方向に直線を伸ばし、各直線がどの計算要素にぶつかるかを判定している。すなわち、ステップＳ３１〜Ｓ３３は、図８のステップＳ２１〜Ｓ２４に相当している。

ステップＳ３３では、計算要素ａからの直線が計算要素Ｂにぶつかるか否かを判定している。この処理の様子を、図１８(ａ)に示す。図１８(ａ)の四角形は、計算要素Ｂとその周囲の計算要素を示す。各四角形内の数値は、各計算要素に立てられたフラグを示す。

フラグ「０」は、初期値に相当する。また、フラグ「１」が立てられた計算要素は、計算要素ａからの直線がその計算要素にぶつかり、かつ、その計算要素がカットオフ角度の範囲内に位置することを意味する。また、フラグ「２」が立てられた計算要素は、計算要素ａからの直線がその計算要素にぶつからないか、その計算要素がカットオフ角度の範囲外に位置することを意味する。

計算要素ａからの直線が計算要素Ｂにぶつかった場合には、計算要素Ｂにフラグ「１」が立てられる（ステップＳ３５）。一方、計算要素ａからの直線が計算要素Ｂにぶつからなかった場合には、計算要素Ｂにフラグ「２」が立てられる（ステップＳ４０）。図８(ａ)は、計算要素Ｂにフラグ「１」が立てられた様子を示している。

なお、ステップＳ３４では、計算要素Ｂが計算要素ａのカットオフ角度の範囲内に位置するか否かを判定している。しかしながら、ステップＳ３１〜Ｓ３２ではカットオフ角度の範囲内の方向のみに直線を伸ばすため、計算要素Ｂは原則的にカットオフ角度の範囲内に位置している。このステップＳ３４は、後に計算要素Ｂの周囲の計算要素についてカットオフ判定を行う際に意味が出てくる。なお、図１７に示すＸ_max、Ｙ_max、Ｘ_min、Ｙ_minはそれぞれ、計算要素Ｂの最大Ｘ座標、最大Ｙ座標、最小Ｘ座標、最小Ｙ座標を示す。

計算要素ａからの直線が計算要素Ｂにぶつかった場合には、図１８(ａ)に示すように、計算要素Ｂに直接的に隣接する各計算要素Ｃについても、計算要素Ｂと同様の処理を行う（ステップＳ３６〜Ｓ３９）。

具体的には、計算要素ａから計算要素Ｃに向けて新たな直線を伸ばし、この直線が計算要素Ｃにぶつかるか否かを判定する（ステップＳ３３）。また、計算要素Ｃが、計算要素ａのカットオフ角度の範囲内に位置するか否かも判定する（ステップＳ３４）。

そして、計算要素ａからの直線がある計算要素Ｃにぶつかり、かつ、その計算要素Ｃがカットオフ角度の範囲内に位置する場合には、その計算要素Ｃにフラグ「１」が立てられる（ステップＳ３５）。一方、計算要素ａからの直線がその計算要素Ｃにぶつからなかった場合や、その計算要素Ｃがカットオフ角度の範囲外に位置する場合には、その計算要素Ｃにフラグ「２」が立てられる（ステップＳ４０）。図８(ｂ)は、各計算要素Ｃにフラグ「１」または「２」が立てられた様子を示している。

ある計算要素Ｃにフラグ「１」が立てられた場合には、図１８(ｃ)に示すように、その計算要素Ｃに直接的に隣接する各計算要素Ｃ’についても、計算要素Ｂと同様の処理を行う（ステップＳ３６〜Ｓ３９）。ただし、すでにフラグ「１」または「２」が立っている計算要素Ｃ’については、この処理は不要である。

本実施形態では、ステップＳ３６〜Ｓ３９の処理を、判定対象の計算要素がなくなるまで繰り返す。具体的には、図１８(ｄ)に示すように、フラグ「１」が立っている計算要素が、フラグ「２」が立っている計算要素に囲まれた状態になるまで、ステップＳ３６〜Ｓ３９の処理が繰り返される。

図１８(ｄ)の符号Ｒは、フラグ「１」が立っている計算要素からなる領域を示す。本実施形態では、この領域Ｒ内に含まれる計算要素を、計算要素ａからの直線がぶつかり、かつ、計算要素ａのカットオフ角度の範囲内に位置する計算要素としてカウントする。

なお、本実施形態では、ステップＳ１２〜Ｓ１４において、各計算要素に固有のローカル座標系を使用したが、代わりに、すべての計算要素に共通のグローバル座標系を使用してもよい。

（４）第１実施形態における計算時間と計算誤差
次に、以上の説明を踏まえて、第１実施形態における計算時間と計算誤差について説明する。

従来の方法では、任意の計算要素Ｂの直接フラックスΓ_B,ne-direct、Γ_B,ion-directを計算するのに、計算要素数Ｎに比例する時間がかかる。理由は、計算要素Ｂに関するループ計算をＮ回繰り返し行うからである。また、従来の方法では、任意の計算要素ａ、Ｂ間の可視係数ν(ａ,Ｂ)や形態係数ｇ(ａ,Ｂ)、ｇ_ionR(ａ,Ｂ)、ｇ_ionS(ａ,Ｂ)を計算するのに、Ｎ²に比例する時間がかかる。理由は、計算要素ａに関するループ計算と計算要素Ｂに関するループ計算を、それぞれＮ回繰り返し行うからである。可視係数や形態係数の計算時間は、鏡面境界条件や周期境界条件を採用するとさらに長くなる。よって、従来の方法における計算時間の多くは、可視係数や形態係数の計算に費やされる。

一方、本実施形態では、図８や図１７に示すように、各計算要素ａから複数の方向に直線を伸ばし、各直線が物質表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定し、この判定結果に基づいて直接フラックス、可視係数、形態係数を計算する。そのため、可視係数と形態係数が、直接フラックスと同様に、計算要素ａに関するループ計算をＮ回繰り返すことで計算される（ステップＳ２２、Ｓ３０を参照）。よって、本実施形態によれば、直接フラックス、可視係数、形態係数の計算時間を、計算要素数Ｎに比例する時間に抑えることができる。

また、本実施形態では、イオン種と中性種の特性の違いを考慮し、反射形態係数ｇ_ionRやスパッタ形態係数ｇ_ionSの計算用にステップＳ２４の判定処理を行う際には、イオン種の反射方向や、イオン種によるスパッタに起因する中性種の発生方向のカットオフ角度を設定し、直線を伸ばす方向をカットオフ角度の範囲内に制限する。さらに、本実施形態では、上記の判定処理を、これらの直線がぶつかった計算要素の周囲の計算要素についても繰り返し適用する。よって、本実施形態によれば、計算の無駄を減らすことで計算時間を短縮することや、無駄な計算の代わりに有用な計算に時間を割くことで計算誤差を低減することが可能となる。なお、中性種の発生量が少ないなど、中性種の発生を無視できる場合には、式（４）におけるスパッタ形態係数ｇ_ionSを含む項を無視して形状シミュレーションを行ってもよい。

また、本実施形態によるｇ、ｇ_ionR、ｇ_ionSの計算によれば、Ｎ²回のループ計算でｇ、ｇ_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）デバイスやディスプレイデバイスなどが挙げられる。

以上、いくつかの実施形態を説明したが、これらの実施形態は、例としてのみ提示したものであり、発明の範囲を限定することを意図したものではない。本明細書で説明した新規な装置、方法、およびプログラムは、その他の様々な形態で実施することができる。また、本明細書で説明した装置、方法、およびプログラムの形態に対し、発明の要旨を逸脱しない範囲内で、種々の省略、置換、変更を行うことができる。添付の特許請求の範囲およびこれに均等な範囲は、発明の範囲や要旨に含まれるこのような形態や変形例を含むように意図されている。

１：シリコン基板、２：シリコン窒化膜、３：シリコン酸化膜、４：貫通孔、
１１：制御部、１２：表示部、１３：入力部、
２１：ＣＰＵ、２２：ＲＯＭ、２３：ＲＡＭ、２４：ＨＤＤ、
２５：メモリドライブ、２６：メモリＩ／Ｆ

Claims (5)

1. 物質の表面を複数の計算要素に分割する分割部と、
   各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定する判定部と、
   前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する計算部とを備え、
   前記判定部は、各計算要素に到達したイオン種が反射する場合における前記形態係数の計算用に前記判定を行う場合、前記イオン種の反射方向のカットオフ角度を設定し、前記直線を伸ばす方向を前記イオン種の反射方向のカットオフ角度の範囲内に制限して前記判定を行い、
   前記判定部は、各計算要素に到達したイオン種によるスパッタにより中性種が発生する場合における前記形態係数の計算用に前記判定を行う場合、前記中性種の発生方向のカットオフ角度を設定し、前記直線を伸ばす方向を前記中性種の発生方向のカットオフ角度の範囲内に制限して前記判定を行い、
   前記判定部は、前記複数の計算要素のうちの第１の計算要素からの直線が第２の計算要素にぶつかった場合には、前記第１の計算要素からの直線が前記第２の計算要素の周囲の第３の計算要素にぶつかるか否かと、前記第３の計算要素が前記第１の計算要素の前記イオン種の反射方向のカットオフ角度の範囲内に位置するか否かとを判定し、
   前記判定部は、前記第３の計算要素として、前記第２の計算要素に隣接する計算要素と、前記第２の計算要素には隣接しないが肯定の判定結果が得られた計算要素には隣接する計算要素とを選択し、選択可能な前記第３の計算要素がなくなるまで前記判定を繰り返し、
   前記判定部は、各計算要素に到達した中性種が各計算要素から再度飛散する場合における前記形態係数の計算用に前記判定を行う場合、前記直線を伸ばす方向のカットオフを行わずに前記判定を行う、
   形状シミュレーション装置。
2. 物質の表面を複数の計算要素に分割する分割部と、
   各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定する判定部と、
   前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する計算部とを備え、
   前記判定部は、各計算要素に到達したイオン種が反射する場合における前記形態係数の計算用に前記判定を行う場合、前記イオン種の反射方向のカットオフ角度を設定し、前記直線を伸ばす方向を前記イオン種の反射方向のカットオフ角度の範囲内に制限して前記判定を行い、
   前記判定部は、前記複数の計算要素のうちの第１の計算要素からの直線が第２の計算要素にぶつかった場合には、前記第１の計算要素からの直線が前記第２の計算要素の周囲の第３の計算要素にぶつかるか否かと、前記第３の計算要素が前記第１の計算要素の前記イオン種の反射方向のカットオフ角度の範囲内に位置するか否かとを判定し、
   前記判定部は、前記第３の計算要素として、前記第２の計算要素に隣接する計算要素と、前記第２の計算要素には隣接しないが肯定の判定結果が得られた計算要素には隣接する計算要素とを選択し、選択可能な前記第３の計算要素がなくなるまで前記判定を繰り返す、
   形状シミュレーション装置。
3. 前記判定部は、各計算要素に到達したイオン種によるスパッタにより中性種が発生する場合における前記形態係数の計算用に前記判定を行う場合、前記中性種の発生方向のカットオフ角度を設定し、前記直線を伸ばす方向を前記中性種の発生方向のカットオフ角度の範囲内に制限して前記判定を行う、請求項２に記載の形状シミュレーション装置。
4. 前記判定部は、各計算要素に到達した中性種が各計算要素から再度飛散する場合における前記形態係数の計算用に前記判定を行う場合、前記直線を伸ばす方向のカットオフを行わずに前記判定を行う、請求項２または３に記載の形状シミュレーション装置。
5. 形状シミュレーション方法をコンピュータに実行させる形状シミュレーションプログラムであって、
   前記形状シミュレーション方法は、
   物質の表面を複数の計算要素に分割し、
   各計算要素から複数の方向に直線を伸ばし、各直線が前記物質の表面にぶつかるか否かと、各直線がどの計算要素にぶつかるかとを判定し、
   前記判定の結果に基づいて、各計算要素に直接的に到達する反応種のフラックスである直接フラックスと、前記計算要素同士の位置関係を示す形態係数とを計算する、
   ことを含み、
   前記判定では、各計算要素に到達したイオン種が反射する場合における前記形態係数の計算用に前記判定を行う場合、前記イオン種の反射方向のカットオフ角度を設定し、前記直線を伸ばす方向を前記イオン種の反射方向のカットオフ角度の範囲内に制限して前記判定を行い、
   前記判定では、前記複数の計算要素のうちの第１の計算要素からの直線が第２の計算要素にぶつかった場合には、前記第１の計算要素からの直線が前記第２の計算要素の周囲の第３の計算要素にぶつかるか否かと、前記第３の計算要素が前記第１の計算要素の前記イオン種の反射方向のカットオフ角度の範囲内に位置するか否かとを判定し、
   前記判定では、前記第３の計算要素として、前記第２の計算要素に隣接する計算要素と、前記第２の計算要素には隣接しないが肯定の判定結果が得られた計算要素には隣接する計算要素とを選択し、選択可能な前記第３の計算要素がなくなるまで前記判定を繰り返す、
   形状シミュレーションプログラム。

Images