JP2734418B2 - デバイス・シミュレーション方法 - Google Patents
デバイス・シミュレーション方法Info
- Publication number
- JP2734418B2 JP2734418B2 JP7208885A JP20888595A JP2734418B2 JP 2734418 B2 JP2734418 B2 JP 2734418B2 JP 7208885 A JP7208885 A JP 7208885A JP 20888595 A JP20888595 A JP 20888595A JP 2734418 B2 JP2734418 B2 JP 2734418B2
- Authority
- JP
- Japan
- Prior art keywords
- electric field
- potential
- lim
- calculation
- mesh
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Classifications
-
- H—ELECTRICITY
- H01—ELECTRIC ELEMENTS
- H01L—SEMICONDUCTOR DEVICES NOT COVERED BY CLASS H10
- H01L22/00—Testing or measuring during manufacture or treatment; Reliability measurements, i.e. testing of parts without further processing to modify the parts as such; Structural arrangements therefor
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- General Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Manufacturing & Machinery (AREA)
- Microelectronics & Electronic Packaging (AREA)
- Power Engineering (AREA)
- Insulated Gate Type Field-Effect Transistor (AREA)
- Testing Or Measuring Of Semiconductors Or The Like (AREA)
- Design And Manufacture Of Integrated Circuits (AREA)
Description
【０００１】
【発明の属する技術分野】本発明はデバイス・シミュレ
ーション方法に関し、特に半導体デバイスの数値解析に
関する。
ーション方法に関し、特に半導体デバイスの数値解析に
関する。
【０００２】
【従来の技術】半導体デバイスの数値解析においては、
キャリア（電子、正孔）を、流体と見做して近似したド
リフト−拡散モデルが広く用いられている。より高次の
近似をしたエネルギー輸送モデルが用いられることもあ
る。本明細書で述べる本発明の手法も、従来手法も、ド
リフト−拡散モデル、エネルギー輸送モデルのどちらに
おいても適用可能である。定常状態でのドリフト−拡散
モデルのデバイス・シミュレーションにおいては、基本
方程式として、以下に示される様な、電荷保存式、電子
電流連続式、正孔電流連続式が設定される。なお、デバ
イス・シミュレーション技術の概説書としては、檀良編
著“プロセス・デバイス・シミュレーション技術”、産
業図書、１９９０年がある。
キャリア（電子、正孔）を、流体と見做して近似したド
リフト−拡散モデルが広く用いられている。より高次の
近似をしたエネルギー輸送モデルが用いられることもあ
る。本明細書で述べる本発明の手法も、従来手法も、ド
リフト−拡散モデル、エネルギー輸送モデルのどちらに
おいても適用可能である。定常状態でのドリフト−拡散
モデルのデバイス・シミュレーションにおいては、基本
方程式として、以下に示される様な、電荷保存式、電子
電流連続式、正孔電流連続式が設定される。なお、デバ
イス・シミュレーション技術の概説書としては、檀良編
著“プロセス・デバイス・シミュレーション技術”、産
業図書、１９９０年がある。
【０００３】
【０００４】μn ：電子移動度 μp ：正孔移動度 Ｄn ：電子拡散係数 Ｄp ：正孔拡散係数 Ｄn ：μn ・ｋB ・Ｔ／ｑ （９） Ｄp ＝μp ・ｋB ・Ｔ／ｑ （１０） ｋB ：ボルツマン定数 Ｔ：温度 上述の式で、解かれるべき変数はポテンシャルΨ、電子
密度ｎ、正孔密度ｐである。一般に、指定された複数の
印加バイアスを境界条件として、順次バイアスを更新し
て、これらの電荷保存式、電子電流連続式、正孔電流連
続式が計算される。
密度ｎ、正孔密度ｐである。一般に、指定された複数の
印加バイアスを境界条件として、順次バイアスを更新し
て、これらの電荷保存式、電子電流連続式、正孔電流連
続式が計算される。
【０００５】これらの電荷保存式、電子電流保存式、正
孔電流電流保存式は非線形な方程式であるため、一般
に、ニュートン法と呼ばれる反復計算を行って、解を求
める。図２０はニュートン法について説明するためのフ
ローチャートである。
孔電流電流保存式は非線形な方程式であるため、一般
に、ニュートン法と呼ばれる反復計算を行って、解を求
める。図２０はニュートン法について説明するためのフ
ローチャートである。
【０００６】変数ｘについて方程式 Ｆ（ｘ）＝０ （１１） が与えられているとする。まず処理手順１０１である初
期値ｘ0 を設定し、処理手順１０２でｉ＝０とおく。あ
る変動量δｘ0 をｘ0 に加えた値が解を与えるならば Ｆ（ｘ0 ＋δｘ0 ）＝０ （１２） である。そして、Ｆ（ｘ）の微係数をＦ′（ｘ0 ）とし
て、Ｆ（ｘ0 ＋δｘ0 ）をδｘ0 について一次のテイラ
ー展開を行うと、 Ｆ（ｘ0 ＋δｘ0 ）≒Ｆ（ｘ0 ）＋Ｆ′（ｘ0 ）δｘ0 ≒０ （１３） そこで処理手順１０３で、 δｘ0 ＝−Ｆ（ｘ0 ）／Ｆ′（ｘ0 ） （１４） 次に、処理手順１０４において、δｘ0 が適当な正の微
小量εより小さいかどうか比較する。もしεより小さけ
れば解が得られたと判断して終了する。もしεより小さ
くなければ処理手順１０５において、 ｘ1 ＝ｘ0 ＋δｘ0 （１５） とおいて、ｘ1 について同様の処理を行なう。一般にｉ
回目の計算でのδｘi がεより小さくなったら（これを
収束したと言い、この判定を収束判定、微小量εを収束
条件という）、そのときのｘi が方程式（１１）の解で
あるとする。ただし、以下の説明では収束したときのｘ
i を最終解と称し、それ以前のｘi を解と称することに
する。
期値ｘ0 を設定し、処理手順１０２でｉ＝０とおく。あ
る変動量δｘ0 をｘ0 に加えた値が解を与えるならば Ｆ（ｘ0 ＋δｘ0 ）＝０ （１２） である。そして、Ｆ（ｘ）の微係数をＦ′（ｘ0 ）とし
て、Ｆ（ｘ0 ＋δｘ0 ）をδｘ0 について一次のテイラ
ー展開を行うと、 Ｆ（ｘ0 ＋δｘ0 ）≒Ｆ（ｘ0 ）＋Ｆ′（ｘ0 ）δｘ0 ≒０ （１３） そこで処理手順１０３で、 δｘ0 ＝−Ｆ（ｘ0 ）／Ｆ′（ｘ0 ） （１４） 次に、処理手順１０４において、δｘ0 が適当な正の微
小量εより小さいかどうか比較する。もしεより小さけ
れば解が得られたと判断して終了する。もしεより小さ
くなければ処理手順１０５において、 ｘ1 ＝ｘ0 ＋δｘ0 （１５） とおいて、ｘ1 について同様の処理を行なう。一般にｉ
回目の計算でのδｘi がεより小さくなったら（これを
収束したと言い、この判定を収束判定、微小量εを収束
条件という）、そのときのｘi が方程式（１１）の解で
あるとする。ただし、以下の説明では収束したときのｘ
i を最終解と称し、それ以前のｘi を解と称することに
する。
【０００７】前述のニュートン法の説明においては、１
個の変数の方程式の場合であったが、デバイス・シミュ
レーションにおいては、例えばＬＤＤ構造ＭＯＳトラン
ジスタのシミュレーションを行なう場合図２１のよう
に、解析領域全体にメッシュを生成し、メッシュ節点上
の変数について方程式を設定する。すなわち、ポテンシ
ャル、電子密度、正孔密度がメッシュ節点の数Ｎだけ変
数として表れるため、３Ｎ個の連立方程式を解くことに
なる。前述の電荷保存式、電子電流連続式、正孔電流連
続式を、右辺の孔を移項した形で、次式のように表す。
個の変数の方程式の場合であったが、デバイス・シミュ
レーションにおいては、例えばＬＤＤ構造ＭＯＳトラン
ジスタのシミュレーションを行なう場合図２１のよう
に、解析領域全体にメッシュを生成し、メッシュ節点上
の変数について方程式を設定する。すなわち、ポテンシ
ャル、電子密度、正孔密度がメッシュ節点の数Ｎだけ変
数として表れるため、３Ｎ個の連立方程式を解くことに
なる。前述の電荷保存式、電子電流連続式、正孔電流連
続式を、右辺の孔を移項した形で、次式のように表す。
【０００８】 ＦΨ （Ψ，ｎ，ｐ）＝０（電荷保存式） （１６） Ｆn （Ψ，ｎ，ｐ）＝０（電子電流連続式） （１７） Ｆp （Ψ，ｎ，ｐ）＝０（正孔電流連続式） （１８） 上式のΨ，ｎ，ｐはそれぞれポテンシャル、電子密度、
正孔密度を表し、またそれぞれＮ個の変数を表す。この
場合、電荷保存式、電子電流連続式、正孔電流連続式を
同時に解くカップル法（結合法）と、電荷保存式、電子
電流連続式、正孔電流連続式を別々に解くガンメル法
（非結合法、またはデカップル法）とがある。
正孔密度を表し、またそれぞれＮ個の変数を表す。この
場合、電荷保存式、電子電流連続式、正孔電流連続式を
同時に解くカップル法（結合法）と、電荷保存式、電子
電流連続式、正孔電流連続式を別々に解くガンメル法
（非結合法、またはデカップル法）とがある。
【０００９】図２２はカップル法について説明するため
のフローチャートである。まず、処理手順１３１で初期
値ｘ0 を設定し、ｉ＝０とおく。ここでｘ0 はΨ0,1 ，
Ψ0,2 ，・・・，Ψ0,N ，ｎ0,1 ，ｎ0,2 ，・・・，ｎ
0,N ，ｐ0,1 ｐ0,2 ，・・・，ｐ0,N なる３Ｎ個の変数
である。添字０は試行回数ｉが０であることに示し、添
字１，２，・・・，ＮはＮ個のメッシュ節点につけた番
号である。
のフローチャートである。まず、処理手順１３１で初期
値ｘ0 を設定し、ｉ＝０とおく。ここでｘ0 はΨ0,1 ，
Ψ0,2 ，・・・，Ψ0,N ，ｎ0,1 ，ｎ0,2 ，・・・，ｎ
0,N ，ｐ0,1 ｐ0,2 ，・・・，ｐ0,N なる３Ｎ個の変数
である。添字０は試行回数ｉが０であることに示し、添
字１，２，・・・，ＮはＮ個のメッシュ節点につけた番
号である。
【００１０】次に、処理手順１３２で Ｆi ′δｘi ＝−Ｆi,i ＝０ （１９） を解く。この式は具体的には次の３Ｎ個の方程式を表わ
す。
す。
【００１１】 Σ（αＦΨi,j ／αΨi,k ）・δΨi,k ＋Σ（αＦΨi,j ／αｎi,k ）・δｎi,k ＋Σ（αＦΨi,j ／αｐi,k ）・δｐi,k ＝−ＦΨi,j （１９−１） Σ（αＦni,j／αΨi,k ）・δΨi,k ＋Σ（αＦni,j／αｎi,k ）・δｎi,k ＋Σ（αＦni,j／αｐi,k ）・δｐi,k ＝−Ｆｎi,j （１９−２） Σ（αＦpi,j／αΨi,k ）・δΨi,k ＋Σ（αＦpi,j／αｎi,k ）・δｎi,k ＋Σ（αＦpi,j／αｐi,k ）・δｐi,k ＝−Ｆｐi,j （１９−３） 式（１９−１）〜（１９−３）において、αは偏微分記
号の代用であり、Σはｋ＝１，２，・・・，Ｎについて
の和をとることを示す。ただし、メッシュ節点ｊ及びこ
れとメッシュエッジを介して隣接しているメッシュ節点
以外のものについては微係数を０とおく。ｊはメッシュ
節点の数Ｎだけある（ｊ＝１，２，３，・・・，Ｎ）。
号の代用であり、Σはｋ＝１，２，・・・，Ｎについて
の和をとることを示す。ただし、メッシュ節点ｊ及びこ
れとメッシュエッジを介して隣接しているメッシュ節点
以外のものについては微係数を０とおく。ｊはメッシュ
節点の数Ｎだけある（ｊ＝１，２，３，・・・，Ｎ）。
【００１２】これらの式よりδΨi,k ，δｎi,k ，δｐ
i,k （ｋ＝１，２，・・・，Ｎ）を得る。
i,k （ｋ＝１，２，・・・，Ｎ）を得る。
【００１３】次に、処理手順１３３で δｘi ＜ε （２０） が成立つかどうか判断する。この式は具体的には次の３
Ｎ個の式を表わす。
Ｎ個の式を表わす。
【００１４】δΨi,k ＜εΨ （２０−１） δｎi,k ＜εn （２０−２） δｐi,k ＜εp （２０−３） εΨ ，εn ，εp は適当な正の微小量である。
【００１５】もしこれが成立すれば最終解が得られたと
して計算を終了（処理手順１３５）する。もし成立しな
ければ、処理手順１３４で ｘi+1 ＝ｘi ＋δｘi （２１） とおく。この式は具体的には次の３Ｎ個の式を表わす。
して計算を終了（処理手順１３５）する。もし成立しな
ければ、処理手順１３４で ｘi+1 ＝ｘi ＋δｘi （２１） とおく。この式は具体的には次の３Ｎ個の式を表わす。
【００１６】 Ψi+1,k ＝Ψi,k ＋δΨi,k （２１−１） ｎi+1,k ＝ｎi,k ＋δｎi,k （２１−２） ｐi+1,k ＝ｐi,k ＋δｐi,k （２１−２） 次にｉ＝ｉ＋１とおいて処理手順１３２以降を繰り返
す。
す。
【００１７】図２３はガンメル法について説明するため
のフローチャートである。処理手順１４１で初期値Ψ
0,j ，ｎ0,j ，ｐ0,j ，ｊ＝１，２，・・・，Ｎを設定
する。
のフローチャートである。処理手順１４１で初期値Ψ
0,j ，ｎ0,j ，ｐ0,j ，ｊ＝１，２，・・・，Ｎを設定
する。
【００１８】次に、処理手順１４２で Ｆ′Ψi δΨi ＝−ＦΨi （２２） を解く。この式は具体的に次の式を表わす。
【００１９】 Σ（αＦΨi,j ／αΨi,k ）・δΨi,k ＝−ＦΨi,j （２３） ここで、αは偏微分記号の代用であり、Σはｋ＝１，
２，・・・，Ｎについての和をとることを示す。ただ
し、メッシュ節点ｊ及びこれとメッシュエッジを介して
隣接しているメッシュ節点以外のものについては微係数
を０とおく。ｊはメッシュ節点の数Ｎだけである。
２，・・・，Ｎについての和をとることを示す。ただ
し、メッシュ節点ｊ及びこれとメッシュエッジを介して
隣接しているメッシュ節点以外のものについては微係数
を０とおく。ｊはメッシュ節点の数Ｎだけである。
【００２０】次に、処理手順１４３で Ｆ′niδｎi ＝−Ｆni （２４） をとく。この式は具体的には次の式を表わす。
【００２１】 Σ（αＦni,j／αｎi,k ）・δｎi,k ＝−Ｆni,j （２５） ここで、αは偏微分記号の代用であり、Σはｋ＝１，
２，・・・，Ｎについての和をとることを示す。メッシ
ュ節点ｊ及びこれとメッシュエッジを介して隣接してい
るメッシュ接点以外のものについては微係数を０とお
く。ｊはメッシュ節点の数だけある。
２，・・・，Ｎについての和をとることを示す。メッシ
ュ節点ｊ及びこれとメッシュエッジを介して隣接してい
るメッシュ接点以外のものについては微係数を０とお
く。ｊはメッシュ節点の数だけある。
【００２２】次に、処理手順１４４で Ｆ′piδｐi ＝−Ｆpi （２６） を解く。この式は具体的には次の式を表わす。
【００２３】 Σ（αＦpi,j／αｐi,k ）・δｐi,k ＝−Ｆpi,j （２７） ここで、αは偏微分記号の代用であり、Σはｋ＝１，
２，・・・，Ｎについての和をとることを示す。メッシ
ュ節点ｊ及びこれとメッシュエッジを介して隣接してい
るメッシュ節点以外のものについては微係数を０とお
く。ｊはメッシュ節点の数だけある。
２，・・・，Ｎについての和をとることを示す。メッシ
ュ節点ｊ及びこれとメッシュエッジを介して隣接してい
るメッシュ節点以外のものについては微係数を０とお
く。ｊはメッシュ節点の数だけある。
【００２４】次に、処理手順１４５で δΨi ＜εΨ （２８） δｎi ＜εn （２９） δｐi ＜εp （３０） が成立つかどうかを判断する。これらの式は具体的には
次の式を表わす。
次の式を表わす。
【００２５】δΨi,k ＜εΨ （３１） δｎi,k ＜εn （３２） δｐi,k ＜εp （３３） ここで、αは偏微分記号の代用であり、εΨ ，εn ，
εp は適当な正の微小量である。
εp は適当な正の微小量である。
【００２６】もし、これらの式が成立てば最終解が得ら
れたとして計算を終了（処理手順１４７）する。もし、
成立たなければ、 Ψi+1 ＝Ψi ＋δΨi （３４） ｎi+1 ＝ｎi ＋δｎi （３５） ｐi+1 ＝ｐi ＋δｐi （３６） とおく。これらの式は具体的には次の３つの式を表わ
す。
れたとして計算を終了（処理手順１４７）する。もし、
成立たなければ、 Ψi+1 ＝Ψi ＋δΨi （３４） ｎi+1 ＝ｎi ＋δｎi （３５） ｐi+1 ＝ｐi ＋δｐi （３６） とおく。これらの式は具体的には次の３つの式を表わ
す。
【００２７】 Ψi+1,k ＝Ψi,k ＋δΨi,k （３７） ｎi+1,k ＝ｎi,k ＋δｎi,k （３８） ｐi+1,k ＝ｐi,k ＋δｐi,k （３９） 次に、ｉ＝ｉ＋１とおいて処理手順１４２以降を繰り返
す。
す。
【００２８】カップル法は、適当な初期値が与えられて
いれば少ない反復回数で最終回を得ることが可能であ
り、初期値依存性が少ないものの反復回数を多く必要と
するガンメル法に比較して、短い計算時間で半導体デバ
イスの解析を行うことが出来る。そのため、カップル法
を用いたデバイスシミュレーションが主流になりつつあ
る。
いれば少ない反復回数で最終回を得ることが可能であ
り、初期値依存性が少ないものの反復回数を多く必要と
するガンメル法に比較して、短い計算時間で半導体デバ
イスの解析を行うことが出来る。そのため、カップル法
を用いたデバイスシミュレーションが主流になりつつあ
る。
【００２９】半導体デバイスのデバイス・シミュレーシ
ョンの基本方程式において、解くべき変数であるポテン
シャル、電子密度、正孔密度の変動領域があまり非線形
性の強くない領域であれば前述のままのニュートン法で
の解の収束が可能であるが、通常はデバイス・シミュレ
ーションの基本方程式は非線形性を有する為、解の変動
量が大きい場合、解が振動して、収束が困難になり、最
終解を得ることが出来なくなる。
ョンの基本方程式において、解くべき変数であるポテン
シャル、電子密度、正孔密度の変動領域があまり非線形
性の強くない領域であれば前述のままのニュートン法で
の解の収束が可能であるが、通常はデバイス・シミュレ
ーションの基本方程式は非線形性を有する為、解の変動
量が大きい場合、解が振動して、収束が困難になり、最
終解を得ることが出来なくなる。
【００３０】前述したような、基本方程式の非線形性に
よる収束の困難を回避するために、解の変動量が所定の
上限値を越えた場合に、これにダンパーを掛けて、反復
を行う手法が一般的に用いられている。
よる収束の困難を回避するために、解の変動量が所定の
上限値を越えた場合に、これにダンパーを掛けて、反復
を行う手法が一般的に用いられている。
【００３１】例えば、ソリッド・ステート・エレクトロ
ニクス（Ｓｏｌｉｄ−ＳｔａｔｅＥｌｅｃｔｒｏｎｉｃ
ｓ）誌、第２７巻、第４号、第３１９頁−第３２８頁、
１９８４年の論文“コンバージェンス・プロパティーズ
・オブ・ニュートンズ・メソッド・フォア・ザ・セミコ
ンダクタ・トランスポート・イウエイションズ・アンド
・ハイブリッド・ソリューション・テクニークス・フォ
ア・マルチディメンショナル・シミュレーション・オブ
・ＶＬＳＩ・デバイス”（Ｃｏｎｖｅｒｇｅｎｃｅ Ｐ
ｒｏｐｅｒｔｉｅｓ ｏｆ Ｎｅｗｔｏｎ’ｓ Ｍｅｔ
ｈｏｄ ｆｏｒ ｔｈｅ Ｓｅｍｉｃｏｎｄｕｃｔｏｒ
Ｔｒａｎｓｐｏｒｔ Ｅｑｕａｔｉｏｎｓ ａｎｄ
Ｈｙｂｒｉｄ Ｓｏｌｕｔｉｏｎ Ｔｅｃｈｎｉｑｕｅ
ｓ ｆｏｒ Ｍｕｌｔｉｄｉｍｅｎｓｉｏｎａｌ Ｓｉ
ｍｕｌａｔｉｏｎ ｏｆ ＶＬＳＩ Ｄｅｖｉｃｅｓ）
にはポテンシャル変動量の上限値を０．７５としてダン
パーを掛けている。このような修正ニュートン法は、そ
の他に、“アイイーイーイー・トランザクションズ・オ
ン・コンピュータ・エイデッド・デザイン”（ＩＥＥＥ
Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｃｏｍｐｕｔａｒ
−Ａｉｄｅｄ Ｄｅｓｉｇｎ）誌，第ＣＡＤ−４巻，第
３号，七月，１９８５年，第１７７頁−第１８９頁や
「アイイーイーイー・トランザクション・オン・エレク
トロン・デバイシス（ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏ
ｎｓ ｏｎ Ｅｌｅｃｔｒｏｎ Ｄｅｖｉｅｓ）誌，第
ＥＤ−３２巻，第１０号，十月、１９８５年、第１９９
２頁−第２００７頁にそれぞれ記載された論文にも使用
されている。
ニクス（Ｓｏｌｉｄ−ＳｔａｔｅＥｌｅｃｔｒｏｎｉｃ
ｓ）誌、第２７巻、第４号、第３１９頁−第３２８頁、
１９８４年の論文“コンバージェンス・プロパティーズ
・オブ・ニュートンズ・メソッド・フォア・ザ・セミコ
ンダクタ・トランスポート・イウエイションズ・アンド
・ハイブリッド・ソリューション・テクニークス・フォ
ア・マルチディメンショナル・シミュレーション・オブ
・ＶＬＳＩ・デバイス”（Ｃｏｎｖｅｒｇｅｎｃｅ Ｐ
ｒｏｐｅｒｔｉｅｓ ｏｆ Ｎｅｗｔｏｎ’ｓ Ｍｅｔ
ｈｏｄ ｆｏｒ ｔｈｅ Ｓｅｍｉｃｏｎｄｕｃｔｏｒ
Ｔｒａｎｓｐｏｒｔ Ｅｑｕａｔｉｏｎｓ ａｎｄ
Ｈｙｂｒｉｄ Ｓｏｌｕｔｉｏｎ Ｔｅｃｈｎｉｑｕｅ
ｓ ｆｏｒ Ｍｕｌｔｉｄｉｍｅｎｓｉｏｎａｌ Ｓｉ
ｍｕｌａｔｉｏｎ ｏｆ ＶＬＳＩ Ｄｅｖｉｃｅｓ）
にはポテンシャル変動量の上限値を０．７５としてダン
パーを掛けている。このような修正ニュートン法は、そ
の他に、“アイイーイーイー・トランザクションズ・オ
ン・コンピュータ・エイデッド・デザイン”（ＩＥＥＥ
Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｃｏｍｐｕｔａｒ
−Ａｉｄｅｄ Ｄｅｓｉｇｎ）誌，第ＣＡＤ−４巻，第
３号，七月，１９８５年，第１７７頁−第１８９頁や
「アイイーイーイー・トランザクション・オン・エレク
トロン・デバイシス（ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏ
ｎｓ ｏｎ Ｅｌｅｃｔｒｏｎ Ｄｅｖｉｅｓ）誌，第
ＥＤ−３２巻，第１０号，十月、１９８５年、第１９９
２頁−第２００７頁にそれぞれ記載された論文にも使用
されている。
【００３２】
【発明が解決しようとする課題】上述した従来の修正ニ
ュートン法では、解の変動量、たとえばポテンシャルの
変動量を反復毎に上限値以下に圧縮して、細かく刻むこ
とによって最終解へ到達出来るようにする。しかし、こ
のような解の変動量の上限値によってダンパーを設定し
た場合でも、基本方程式の非線形性は一様でない為、非
線形性の強くない領域での解の振動に合わせて設定した
上限値では、非線形性の強くなる領域において収束が困
難になる。逆に、強い非線形性領域に合わせて上限値を
設定したのでは、非線形性の強くない反復回数を不必要
に増大させ、計算時間が長くなってしまう。
ュートン法では、解の変動量、たとえばポテンシャルの
変動量を反復毎に上限値以下に圧縮して、細かく刻むこ
とによって最終解へ到達出来るようにする。しかし、こ
のような解の変動量の上限値によってダンパーを設定し
た場合でも、基本方程式の非線形性は一様でない為、非
線形性の強くない領域での解の振動に合わせて設定した
上限値では、非線形性の強くなる領域において収束が困
難になる。逆に、強い非線形性領域に合わせて上限値を
設定したのでは、非線形性の強くない反復回数を不必要
に増大させ、計算時間が長くなってしまう。
【００３３】従って本発明の目的は、解の収束を確保で
きかつ計算時間の短縮が可能なデバイス・シミュレーシ
ョン法を提供することにある。
きかつ計算時間の短縮が可能なデバイス・シミュレーシ
ョン法を提供することにある。
【００３４】
【課題を解決するための手段】本発明第１のデバイス・
シミュレーション方法は、半導体デバイスの数値解析の
ニュートン反復法による計算の過程で、ｉ回目の計算結
果から得られるポテンシャルΨi+1 を用いて電界強度Ｅ
i+1 を計算し、前記電界強度の絶対値｜Ｅi+1｜が所定
の上限値Ｅlim を越えた場合、更に電界変動量δＥi ＝
Ｅi+1 −Ｅi の絶対値｜δＥi ｜が所定の上限値δＥ
lim を越えたときに、ポテンシャル変動量δΨi ＝Ψ
i+1 −Ψi に所定のダンパーＤ1iを掛けたＤ1i・δΨi
とΨi との和Ψi+1 a ＝Ψi ＋Ｄ1i・δΨi をΨi+1 の
代りに使用して（ｉ＋１）回目の計算を行なう手順を有
するというものである。
シミュレーション方法は、半導体デバイスの数値解析の
ニュートン反復法による計算の過程で、ｉ回目の計算結
果から得られるポテンシャルΨi+1 を用いて電界強度Ｅ
i+1 を計算し、前記電界強度の絶対値｜Ｅi+1｜が所定
の上限値Ｅlim を越えた場合、更に電界変動量δＥi ＝
Ｅi+1 −Ｅi の絶対値｜δＥi ｜が所定の上限値δＥ
lim を越えたときに、ポテンシャル変動量δΨi ＝Ψ
i+1 −Ψi に所定のダンパーＤ1iを掛けたＤ1i・δΨi
とΨi との和Ψi+1 a ＝Ψi ＋Ｄ1i・δΨi をΨi+1 の
代りに使用して（ｉ＋１）回目の計算を行なう手順を有
するというものである。
【００３５】また、本発明の第２のデバイス・シミュレ
ーション方法は、半導体デバイスの数値解析のニュート
ン反復法による計算の過程で、ｉ回目の計算結果のポテ
ンシャルΨi を用いて電界強度Ｅi を計算し、前記電界
強度の絶対値｜Ｅi ｜が所定の上限値Ｅlim を越えた場
合、更にｉ回目の計算結果から得られるポテンシャルΨ
i+1 を用いて電界強度Ｅi+1 を計算し、電界強度比の絶
対値｜Ｅi+1 ／Ｅi ｜が所定値Ｒlim を越えたときに、
ポテンシャル変動量δΨi ＝Ψi+1 −Ψi に所定のダン
パーＤ2iをかけたＤ2i・δΨi とΨi との和Ψi+1 b ＝
Ψi ＋Ｄ2i・δΨi をΨi+1 の代りに使用して（ｉ＋
１）回目の計算を行なう手順を有する。これらの場合、
電荷保存式、電子電流連続式及び正孔電流連続式を所定
の境界条件のもとで好ましくはカップル法で解くことが
できる。
ーション方法は、半導体デバイスの数値解析のニュート
ン反復法による計算の過程で、ｉ回目の計算結果のポテ
ンシャルΨi を用いて電界強度Ｅi を計算し、前記電界
強度の絶対値｜Ｅi ｜が所定の上限値Ｅlim を越えた場
合、更にｉ回目の計算結果から得られるポテンシャルΨ
i+1 を用いて電界強度Ｅi+1 を計算し、電界強度比の絶
対値｜Ｅi+1 ／Ｅi ｜が所定値Ｒlim を越えたときに、
ポテンシャル変動量δΨi ＝Ψi+1 −Ψi に所定のダン
パーＤ2iをかけたＤ2i・δΨi とΨi との和Ψi+1 b ＝
Ψi ＋Ｄ2i・δΨi をΨi+1 の代りに使用して（ｉ＋
１）回目の計算を行なう手順を有する。これらの場合、
電荷保存式、電子電流連続式及び正孔電流連続式を所定
の境界条件のもとで好ましくはカップル法で解くことが
できる。
【００３６】解が振動して収束しない理由は、高電界領
域でこれらの方程式の非線形性が強いことに基いてい
る。従って、このような高電界領域でのポテンシャル変
動量にダンパーを掛けることによって収束させることが
できる。高電界領域以外ではダンパーを掛けないので反
復回数の増大を回避できる。
域でこれらの方程式の非線形性が強いことに基いてい
る。従って、このような高電界領域でのポテンシャル変
動量にダンパーを掛けることによって収束させることが
できる。高電界領域以外ではダンパーを掛けないので反
復回数の増大を回避できる。
【００３７】
【発明の実施の形態】図１，図２は本発明の第１の実施
の形態について説明するためのフローチャートの部分図
である。
の形態について説明するためのフローチャートの部分図
である。
【００３８】方程式（１６），（１７）及び（１８）を
所定の境界条件のものでカップル法により解く。図１に
おいて、処理手順１，２，３は図２２の従来の処理手順
１３１，１３２，１３３と全く同様である。処理手順３
で式（２０−１），（２０−２），（２０−３）が成立
すれば図２の端子Ｄを経て終了（処理手順１２）する。
もし、成立しなければ、処理手順４を行ない、 ｘi+1 ＝ｘi ＋δｘi （２１） とおく。この式は具体的には次の３Ｎ個の式を表わす。
所定の境界条件のものでカップル法により解く。図１に
おいて、処理手順１，２，３は図２２の従来の処理手順
１３１，１３２，１３３と全く同様である。処理手順３
で式（２０−１），（２０−２），（２０−３）が成立
すれば図２の端子Ｄを経て終了（処理手順１２）する。
もし、成立しなければ、処理手順４を行ない、 ｘi+1 ＝ｘi ＋δｘi （２１） とおく。この式は具体的には次の３Ｎ個の式を表わす。
【００３９】 Ψi+1,k ＝Ψi,k ＋δΨi,k （２１−１） ｎi+1,k ＝ｎi,k ＋δｎi,k （２１−２） ｐi+1,k ＝ｐi,k ＋δｐi,k （２１−３） 次に、処理手順５で電界強度Ｅi+1 を計算する。すなわ
ち、 Ｅi+1,l ＝（Ψi+1,j −Ψi+1,k ）／ｄl （４０） ポテンシャルΨi+1 はメッシュ節点上で定義されるが、
電界強度Ｅi+1 は基本方程式から直接得られる物理量で
はなく、ポテンシャルΨi+1 を用いて計算される。メッ
シュエッジをｌで表わし、その長さをｄｌとし、その両
端のメッシュ節点をｊ，ｋで表わすと、メッシュエッジ
上の電界強度は式（４０）で定義される。
ち、 Ｅi+1,l ＝（Ψi+1,j −Ψi+1,k ）／ｄl （４０） ポテンシャルΨi+1 はメッシュ節点上で定義されるが、
電界強度Ｅi+1 は基本方程式から直接得られる物理量で
はなく、ポテンシャルΨi+1 を用いて計算される。メッ
シュエッジをｌで表わし、その長さをｄｌとし、その両
端のメッシュ節点をｊ，ｋで表わすと、メッシュエッジ
上の電界強度は式（４０）で定義される。
【００４０】次に、処理手順６において、電界強度Ｅ
i+1 の絶対値があらかじめ適当に設定された上限値Ｅ
lim より大きくなるかどうかを判断する。すなわち、 ｜Ｅi+1,l ｜＞Ｅlim （４１） を満たすメッシュエッジｌが存在する（高電界領域）か
どうかを判断する。もし、このようなメッシュエッジｌ
がなければ、端子Ｂを経て図２の処理手順１１に進みｉ
＝ｉ＋１とおき、端子Ａを経て処理手順２を行なう。も
し、このようなメッシュエッジｌが１つでもあれば（高
電界領域ならば）端子Ｃを経て図２の処理手順７に進
み、前回反復で得られたポテンシャルΨｉを用いて、 Ｅi,l ＝（Ψi,j −Ψi,k ）／ｄl （４２） により電界強度Ｅi を求める。次に、処理手順８で、 δＥi,l ＝Ｅi+1,l −Ｅi,l （４３） により、電界変動量δＥi を求める。次に、処理手順９
で、適当に設定された電界強度の上限値δＥlim よりも
δＥi の絶対値が大きいかどうかを判断する。すなわち ｜δＥi,l ｜＞δＥlim （４４） なるメッシュエッジｌがあるかどうかを判断する。も
し、このようなメッシュエジｌが存在しないならば処理
手順１１に進む。もし、存在するならば処理手順１０に
進み、ポテンシャル変動量δΨi ＝Ψi+1 −Ψi にダン
パーＤ1iをかける。すなわち、 δΨi,l a ＝Ｄ1i・δΨi,l （４５） Ｄ1i＝δＥlim ／｜Ｅi+1,M −Ｅi,M ｜ （４６） で定義されるδΨi,l a をδΨi,l の代りに使用して、
処理時１１に進む。ただし、Ｅi+1,M ，Ｅi,M はδＥ
i,l の絶対値｜δＥi,l ｜が最大値｜δＥi,MAX ｜とな
るメッシュエッジＭでの電界強度を表わす。従って、メ
ッシュエッジＭの長さをｄM 、その両端のメッシュ節点
をｊ，ｋで表わすと、 δＥi,M a ＝（δΨi,j a −δΨi,k a ）／ｄM （４７） ＝Ｄ1i・δＥi,M （４８） となる。修正された電界変動量δＥi,l a をＤ1i・δＥ
i,l で定義すると、 すなわち、式（４５），（４６）で定義される修正ポテ
ンシャル変動量δΨi,l a をδΨi,l の代りに使用する
と、あたかも電界変動量の絶対値が上限値以下の領域で
計算を行なうことができる。
i+1 の絶対値があらかじめ適当に設定された上限値Ｅ
lim より大きくなるかどうかを判断する。すなわち、 ｜Ｅi+1,l ｜＞Ｅlim （４１） を満たすメッシュエッジｌが存在する（高電界領域）か
どうかを判断する。もし、このようなメッシュエッジｌ
がなければ、端子Ｂを経て図２の処理手順１１に進みｉ
＝ｉ＋１とおき、端子Ａを経て処理手順２を行なう。も
し、このようなメッシュエッジｌが１つでもあれば（高
電界領域ならば）端子Ｃを経て図２の処理手順７に進
み、前回反復で得られたポテンシャルΨｉを用いて、 Ｅi,l ＝（Ψi,j −Ψi,k ）／ｄl （４２） により電界強度Ｅi を求める。次に、処理手順８で、 δＥi,l ＝Ｅi+1,l −Ｅi,l （４３） により、電界変動量δＥi を求める。次に、処理手順９
で、適当に設定された電界強度の上限値δＥlim よりも
δＥi の絶対値が大きいかどうかを判断する。すなわち ｜δＥi,l ｜＞δＥlim （４４） なるメッシュエッジｌがあるかどうかを判断する。も
し、このようなメッシュエジｌが存在しないならば処理
手順１１に進む。もし、存在するならば処理手順１０に
進み、ポテンシャル変動量δΨi ＝Ψi+1 −Ψi にダン
パーＤ1iをかける。すなわち、 δΨi,l a ＝Ｄ1i・δΨi,l （４５） Ｄ1i＝δＥlim ／｜Ｅi+1,M −Ｅi,M ｜ （４６） で定義されるδΨi,l a をδΨi,l の代りに使用して、
処理時１１に進む。ただし、Ｅi+1,M ，Ｅi,M はδＥ
i,l の絶対値｜δＥi,l ｜が最大値｜δＥi,MAX ｜とな
るメッシュエッジＭでの電界強度を表わす。従って、メ
ッシュエッジＭの長さをｄM 、その両端のメッシュ節点
をｊ，ｋで表わすと、 δＥi,M a ＝（δΨi,j a −δΨi,k a ）／ｄM （４７） ＝Ｄ1i・δＥi,M （４８） となる。修正された電界変動量δＥi,l a をＤ1i・δＥ
i,l で定義すると、 すなわち、式（４５），（４６）で定義される修正ポテ
ンシャル変動量δΨi,l a をδΨi,l の代りに使用する
と、あたかも電界変動量の絶対値が上限値以下の領域で
計算を行なうことができる。
【００４１】前述したように、ニュートン法で基本方程
式を解く場合に振動して収束が困難になる原因はデバイ
ス・シミュレーションの基本方程式の非線形性にある。
基本方程式の非線形性は移動度の電界強度依存性による
ところが大きい。例えば、その修正モデルが広く用いら
れているシャーフェッタ・ガンメル（Ｓｃｈａｒｆｅｔ
ｔｅｒ ａｎｄ Ｇｕｍｍｅｌ）の移動度モデルを図３
に示す。これは、正孔移動度について、横軸に電界強度
を取って、電界強度依存性の様子を表している。電界が
１．０×１０4 Ｖ／ｃｍを越えるような高電界領域で
は、移動度が急激に低下していて、非線形性が強い。そ
のため、このような非線形性の強い高電界領域では、わ
ずかな電界強度の変動によって、移動度が大きく変化す
る。その結果、デバイス・シミュレーションの基本方程
式を解く上では、電流連続式の中に移動度が取り込まれ
ているために、このような非線形性の強い高電界領域で
は、移動度の大きな変動によって、解が振動してしま
い、収束しない場合がある。
式を解く場合に振動して収束が困難になる原因はデバイ
ス・シミュレーションの基本方程式の非線形性にある。
基本方程式の非線形性は移動度の電界強度依存性による
ところが大きい。例えば、その修正モデルが広く用いら
れているシャーフェッタ・ガンメル（Ｓｃｈａｒｆｅｔ
ｔｅｒ ａｎｄ Ｇｕｍｍｅｌ）の移動度モデルを図３
に示す。これは、正孔移動度について、横軸に電界強度
を取って、電界強度依存性の様子を表している。電界が
１．０×１０4 Ｖ／ｃｍを越えるような高電界領域で
は、移動度が急激に低下していて、非線形性が強い。そ
のため、このような非線形性の強い高電界領域では、わ
ずかな電界強度の変動によって、移動度が大きく変化す
る。その結果、デバイス・シミュレーションの基本方程
式を解く上では、電流連続式の中に移動度が取り込まれ
ているために、このような非線形性の強い高電界領域で
は、移動度の大きな変動によって、解が振動してしま
い、収束しない場合がある。
【００４２】本発明では、高電界領域でポテンシャル変
動量にダンパーをかけてニュートン法を修正することに
よって解の収束を確保し、高電界領域以外では通常のニ
ュートン法を使用することによって非線形性の強くない
領域での不必要な反復回数の増大を回避できる。
動量にダンパーをかけてニュートン法を修正することに
よって解の収束を確保し、高電界領域以外では通常のニ
ュートン法を使用することによって非線形性の強くない
領域での不必要な反復回数の増大を回避できる。
【００４３】図４，図５は本発明の第２の実施の形態に
ついて説明するためのフローチャートの部分図である。
ついて説明するためのフローチャートの部分図である。
【００４４】図４の処理手順１３，１４，１５は図１の
処理手順１，２，３と同じである。処理手順１５で式
（２０−１），（２０−２），（２０−３）が成立すれ
ば図５の端子Ｄを経て終了（処理手順２４）する。も
し、成立しなければ処理手順１６でポテンシャルΨ
i （今回反復の処理手順１４で求めた）から電界強度Ｅ
i を計算する。
処理手順１，２，３と同じである。処理手順１５で式
（２０−１），（２０−２），（２０−３）が成立すれ
ば図５の端子Ｄを経て終了（処理手順２４）する。も
し、成立しなければ処理手順１６でポテンシャルΨ
i （今回反復の処理手順１４で求めた）から電界強度Ｅ
i を計算する。
【００４５】Ψi ＝Ψi-1 ＋δΨi-1 （５１） これは、具体的には次のＮ個の式を表わす。
【００４６】 Ψi,j ＝Ψi-1,j ＋δΨi-1,j （５２） ｊはメッシュ節点の数Ｎだけある。電界強度は基本方程
式から直接得られる物理量ではなく、ポテンシャルから
得られる量であり、メッシュエッジ上で定義される。す
なわち、メッシュエッジをｌで表わし、その長さをｄl
とし、その両端のメッシュ節点をｊ，ｋで表すと、メッ
シュエッジ上の電界強度Ｅi,l は Ｅi,l ＝（Ψi,j −Ψi,k ）／ｄl （５３） で定義される。次に処理手順１７で、電界強度Ｅi の絶
対値があらかじめ適当に設定された上限値Ｅlim より大
きくなるかどうかを判断する。すなわち、 ｜Ｅi,l ｜＞Ｅlim （５４） を満たすメッシュエッジｌがあるかどうかを判断する。
もし、このようなメッシュエッジがなければ図５の端子
Ｂを経て処理手順２３に進む。もしあれば、処理手順１
８へ進み ｘi+1 ＝ｘi ＋δｘi （５５） とおく。この式は具体的には次の３Ｎ個の式を表わす。
式から直接得られる物理量ではなく、ポテンシャルから
得られる量であり、メッシュエッジ上で定義される。す
なわち、メッシュエッジをｌで表わし、その長さをｄl
とし、その両端のメッシュ節点をｊ，ｋで表すと、メッ
シュエッジ上の電界強度Ｅi,l は Ｅi,l ＝（Ψi,j −Ψi,k ）／ｄl （５３） で定義される。次に処理手順１７で、電界強度Ｅi の絶
対値があらかじめ適当に設定された上限値Ｅlim より大
きくなるかどうかを判断する。すなわち、 ｜Ｅi,l ｜＞Ｅlim （５４） を満たすメッシュエッジｌがあるかどうかを判断する。
もし、このようなメッシュエッジがなければ図５の端子
Ｂを経て処理手順２３に進む。もしあれば、処理手順１
８へ進み ｘi+1 ＝ｘi ＋δｘi （５５） とおく。この式は具体的には次の３Ｎ個の式を表わす。
【００４７】 Ψi+1,k ＝Ψi,k ＋δΨi,k （５５−１） ｎi+1,k ＝ｎi,k ＋δｎi,k （５５−２） ｐi+1,k ＝ｐi,k ＋δｐi,k （５５−３） 次に、処理手順１９で電界強度Ｅi+1 を計算する。すな
わち、 Ｅi+1,l ＝（Ψi+1,j −Ψi+1,k ）／ｄl （５６） 次に図５の節点Ｃを経て処理手順２０で電界強度比の絶
対値｜Ｅi+1 ／Ｅi ｜を計算する。すなわち、 Ｒi,l ＝｜Ｅi+1,l ／Ｅi,l ｜ （５７） 分母にＥi を用いることによって自動的にゼロになる割
り算を防ぐことができる。
わち、 Ｅi+1,l ＝（Ψi+1,j −Ψi+1,k ）／ｄl （５６） 次に図５の節点Ｃを経て処理手順２０で電界強度比の絶
対値｜Ｅi+1 ／Ｅi ｜を計算する。すなわち、 Ｒi,l ＝｜Ｅi+1,l ／Ｅi,l ｜ （５７） 分母にＥi を用いることによって自動的にゼロになる割
り算を防ぐことができる。
【００４８】次に、処理手順２１で、このＲi,l の絶対
値が適当に設定された上限値Ｒlimより大きいかどうか
を判断する。すなわち、 Ｒi,l ＞Ｒlim （５４） が成立つメッシュエッジｌがあるかどうかを判断する。
もしなければ処理手順２３へ進む。もしあれば、処理手
順２２でポテンシャル変動量δΨi にダンパーＤ2iを掛
ける。すなわち、 δΨi,l b ＝Ｄ2i・δΨi,l （５５） Ｄ2i＝｜Ｅi,M ／δＥi,M ｜・（Ｒlim −１） （５６） で定義されるδΨi,l b をδΨi,l の代りに使用して処
理手順２３に進む。ただし、Ｅi,M は電界強度比の絶対
値が最大値Ｒi,MAX となるメッシュエッジＭ上での電界
強度を表わす。従って、メッシュエッジＭの長さを
ｄM 、その両端のメッシュ節点をｊ，ｋで表わすと δＥi,M b ＝（δΨi,j b −δΨi,k b ）／ｄM （５７） ＝Ｄ2i・δＥi,M （５８） となる。そこで、Ｒi,MAX b を次式で定義する。
値が適当に設定された上限値Ｒlimより大きいかどうか
を判断する。すなわち、 Ｒi,l ＞Ｒlim （５４） が成立つメッシュエッジｌがあるかどうかを判断する。
もしなければ処理手順２３へ進む。もしあれば、処理手
順２２でポテンシャル変動量δΨi にダンパーＤ2iを掛
ける。すなわち、 δΨi,l b ＝Ｄ2i・δΨi,l （５５） Ｄ2i＝｜Ｅi,M ／δＥi,M ｜・（Ｒlim −１） （５６） で定義されるδΨi,l b をδΨi,l の代りに使用して処
理手順２３に進む。ただし、Ｅi,M は電界強度比の絶対
値が最大値Ｒi,MAX となるメッシュエッジＭ上での電界
強度を表わす。従って、メッシュエッジＭの長さを
ｄM 、その両端のメッシュ節点をｊ，ｋで表わすと δＥi,M b ＝（δΨi,j b −δΨi,k b ）／ｄM （５７） ＝Ｄ2i・δＥi,M （５８） となる。そこで、Ｒi,MAX b を次式で定義する。
【００４９】 Ｒi,MAX b ＝｜Ｅi+1,M b ／Ｅi,M ｜ （５９） ＝｜（Ｅi,M ＋Ｄ2i・δＥi,M ）／Ｅi,M ｜ （６０） ＝Ｒlim （６１） すなわち、式（５５），（５６）で定義される修正ポテ
ンシャル変動量δΨi,l b をδΨi,l を代りに使用する
と、あたかも電界強度比の絶対値が上限値以下の領域で
計算を行なうことができる。
ンシャル変動量δΨi,l b をδΨi,l を代りに使用する
と、あたかも電界強度比の絶対値が上限値以下の領域で
計算を行なうことができる。
【００５０】第１の実施の形態と同様に解の収束の確保
と計算時間の不必要な増大を回避できるが、第２の実施
の形態では高電界領域の広い領域で適切なダンパーを設
定できる点で優れている。第１の実施の形態で、例え
ば、δＥlim が２．０×１０4Ｖ／ｃｍに設定されてい
るとき、Ｅi,M が１．０×１０4 Ｖ／ｃｍ程度の領域で
は適切かもしれないが、Ｅi,M が１．０×１０5 Ｖ／ｃ
ｍ程度の領域では上限値としては小さすぎて、反復回数
を増大させる恐れがある。一方、第２の実施の形態なら
ば、電界変動比で上限値を見ているので、この場合にも
適切にダンパーを掛けることが可能である。
と計算時間の不必要な増大を回避できるが、第２の実施
の形態では高電界領域の広い領域で適切なダンパーを設
定できる点で優れている。第１の実施の形態で、例え
ば、δＥlim が２．０×１０4Ｖ／ｃｍに設定されてい
るとき、Ｅi,M が１．０×１０4 Ｖ／ｃｍ程度の領域で
は適切かもしれないが、Ｅi,M が１．０×１０5 Ｖ／ｃ
ｍ程度の領域では上限値としては小さすぎて、反復回数
を増大させる恐れがある。一方、第２の実施の形態なら
ば、電界変動比で上限値を見ているので、この場合にも
適切にダンパーを掛けることが可能である。
【００５１】
【実施例】次に、本発明の実施例について説明する。
【００５２】図６は第１の実施例について説明するため
のＭＯＳトランジスタを示す断面図で、このＭＯＳトラ
ンジスタは、Ｐ型シリコン基板３１の表面にゲート酸化
膜３５を介してゲート電極Ｇが設けられている。３２は
ソース領域、３３はドレイン領域でいずれもＮ型拡散層
で、破線はＰＮ接合位置を示す。この破線はプロセス・
シミュレーションの結果をそのまま図示してある。３４
はチャネルドープ領域、３６は酸化シリコン膜でなる側
壁スペーサ、Ｓはソース電極、Ｄはドレイン電極であ
る。図７はこのようなＭＯＳトランジスタのゲート酸化
膜界面における不純物濃度分布を示し、実線はドナー
の、破線はアクセプタの分布をそれぞれ示す。図８に生
成したメッシュを示す。黒く塗りつぶされたように見え
るところはメッシュが細いところである。図９は本実施
例で使用した電子移動度モデル、図１０は正孔移動度モ
デルである。
のＭＯＳトランジスタを示す断面図で、このＭＯＳトラ
ンジスタは、Ｐ型シリコン基板３１の表面にゲート酸化
膜３５を介してゲート電極Ｇが設けられている。３２は
ソース領域、３３はドレイン領域でいずれもＮ型拡散層
で、破線はＰＮ接合位置を示す。この破線はプロセス・
シミュレーションの結果をそのまま図示してある。３４
はチャネルドープ領域、３６は酸化シリコン膜でなる側
壁スペーサ、Ｓはソース電極、Ｄはドレイン電極であ
る。図７はこのようなＭＯＳトランジスタのゲート酸化
膜界面における不純物濃度分布を示し、実線はドナー
の、破線はアクセプタの分布をそれぞれ示す。図８に生
成したメッシュを示す。黒く塗りつぶされたように見え
るところはメッシュが細いところである。図９は本実施
例で使用した電子移動度モデル、図１０は正孔移動度モ
デルである。
【００５３】本発明の第２の実施の形態の手法を用いて
デバイス・シミュレーションを行なった。
デバイス・シミュレーションを行なった。
【００５４】図１１は、ソース電極Ｓに０．０Ｖ、ゲー
ト電極Ｇに１．５Ｖ、ドレイン電極Ｄに０．５Ｖ、基板
電極Ｂに０．０Ｖのバイアスを印加した状態を解析した
後、ゲート電極印加バイアスのみを２．０Ｖにした状態
を解析しようとした時の、収束の様子を示したものであ
る。横軸に反復回数をとり、横軸にポテンサル変動量の
最大値の絶対値｜δΨMAX ｜をとっている。
ト電極Ｇに１．５Ｖ、ドレイン電極Ｄに０．５Ｖ、基板
電極Ｂに０．０Ｖのバイアスを印加した状態を解析した
後、ゲート電極印加バイアスのみを２．０Ｖにした状態
を解析しようとした時の、収束の様子を示したものであ
る。横軸に反復回数をとり、横軸にポテンサル変動量の
最大値の絶対値｜δΨMAX ｜をとっている。
【００５５】従来手法ではポテンシャル変動量の最大値
の絶対置｜δΨMAX ｜が減少せず、収束が困難であるこ
とがわかる。ここでは、従来手法でのポテンシャル変動
量の上限値は、０．０７５Ｖとした。一方、本発明の手
法ではすみやかにポテンシャル変動量が減少し、解の収
束を得ていることがわかる。ここでは、電界変動量比の
上限値Ｒlim を３とし、Ｅlim は１．０×１０4 Ｖ／ｃ
ｍとした。
の絶対置｜δΨMAX ｜が減少せず、収束が困難であるこ
とがわかる。ここでは、従来手法でのポテンシャル変動
量の上限値は、０．０７５Ｖとした。一方、本発明の手
法ではすみやかにポテンシャル変動量が減少し、解の収
束を得ていることがわかる。ここでは、電界変動量比の
上限値Ｒlim を３とし、Ｅlim は１．０×１０4 Ｖ／ｃ
ｍとした。
【００５６】図１２は、本実施例の結果得られたシリコ
ン基板表面のポテンシャル分布を示し、図１３は同じく
電子移動度の分布を示す。
ン基板表面のポテンシャル分布を示し、図１３は同じく
電子移動度の分布を示す。
【００５７】次に、第２の実施例について説明する。
【００５８】図１４は第２の実施例について説明するた
めのＭＯＳトランジスタを示す断面図、図１５はゲート
酸化膜界面における不純物濃度分布を示し、実線はドナ
ーの、破線はアクセプタの分布をそれぞれ示す。図１６
は生成したメッシュを示し、黒く塗りつぶしたように見
えるところはメッシュが細いところである。電子移動度
モデル、正孔移動度モデルは第１の実施例と同じものを
使用した。
めのＭＯＳトランジスタを示す断面図、図１５はゲート
酸化膜界面における不純物濃度分布を示し、実線はドナ
ーの、破線はアクセプタの分布をそれぞれ示す。図１６
は生成したメッシュを示し、黒く塗りつぶしたように見
えるところはメッシュが細いところである。電子移動度
モデル、正孔移動度モデルは第１の実施例と同じものを
使用した。
【００５９】本発明の第２の実施の形態の手法を用いて
デバイス・シミュレーションを行なった。ただし、式
（５），（６）のキャリア生成項に局所的インパクトイ
オン化モデルを用いた。すなわち、
デバイス・シミュレーションを行なった。ただし、式
（５），（６）のキャリア生成項に局所的インパクトイ
オン化モデルを用いた。すなわち、
【００６０】
【００６１】なお、インパクトイオン化モデルは「電子
通信学会技術研究報告」誌、第８４巻、第１８０号、第
４５頁−第５２頁、１９８４年に紹介されている。
通信学会技術研究報告」誌、第８４巻、第１８０号、第
４５頁−第５２頁、１９８４年に紹介されている。
【００６２】図１７は、ソース電極Ｓに０．０Ｖ、ゲー
ト電極Ｇに０．５Ｖ、ドレイン電極Ｄに０．０Ｖ、基板
電極Ｂに０．０Ｖのバイアスを印加した状態を解析した
後、ゲート電極印加バイアスのみを１．０Ｖにした状態
を解析しようとした時の収束の様子を示したものであ
る。横軸に反復回数をとり、縦軸にポテンシャル変動量
の最大値の絶対値｜δΨMAX ｜をとっている。
ト電極Ｇに０．５Ｖ、ドレイン電極Ｄに０．０Ｖ、基板
電極Ｂに０．０Ｖのバイアスを印加した状態を解析した
後、ゲート電極印加バイアスのみを１．０Ｖにした状態
を解析しようとした時の収束の様子を示したものであ
る。横軸に反復回数をとり、縦軸にポテンシャル変動量
の最大値の絶対値｜δΨMAX ｜をとっている。
【００６３】従来手法ではポテンシャル変動量の増大値
の絶対値｜δΨMAX ｜が減少せず、収束が困難であるこ
とがわかる。ここでは、従来手法でのポテンシャル変動
量の上限値は、０．０７５Ｖとした。一方、本発明の手
法ではすみやかにポテンシャル変動量が減少し、解の収
束を得ていることがわかる。ここでは、電界変動比の上
限値Ｒlim を３とし、Ｅlim は１．０×１０4 Ｖ／ｃｍ
とした。
の絶対値｜δΨMAX ｜が減少せず、収束が困難であるこ
とがわかる。ここでは、従来手法でのポテンシャル変動
量の上限値は、０．０７５Ｖとした。一方、本発明の手
法ではすみやかにポテンシャル変動量が減少し、解の収
束を得ていることがわかる。ここでは、電界変動比の上
限値Ｒlim を３とし、Ｅlim は１．０×１０4 Ｖ／ｃｍ
とした。
【００６４】図１８は本実施例の結果得られたシリコン
基板表面のポテンシャル分布を示し、図１９は同じくイ
ンパクトイオン化率の分布を示す。
基板表面のポテンシャル分布を示し、図１９は同じくイ
ンパクトイオン化率の分布を示す。
【００６５】以上、半導体デバイスの数値解析のニュー
トン反復法の効率的な修正法についてカップル法を例に
あげて説明したが、本発明はガンメル法にも適用可能で
あることは当業者にとって明らかであろう。また、ドリ
フト−拡散モデルの場合について説明したが、エネルギ
ー輸送モデルの場合にも本発明を適用することができ
る。
トン反復法の効率的な修正法についてカップル法を例に
あげて説明したが、本発明はガンメル法にも適用可能で
あることは当業者にとって明らかであろう。また、ドリ
フト−拡散モデルの場合について説明したが、エネルギ
ー輸送モデルの場合にも本発明を適用することができ
る。
【００６６】数値解析の結果は、計算に用いたコンピュ
ータの記憶装置を経てハードディスクなどのファイルに
格納される。このハードディスクの内容はしかるべきプ
ログラムを使用して図１２，図１３などのようなグラフ
としてディスプレイに表示させたり、プリントアウトす
ることができることはいうまでもない。
ータの記憶装置を経てハードディスクなどのファイルに
格納される。このハードディスクの内容はしかるべきプ
ログラムを使用して図１２，図１３などのようなグラフ
としてディスプレイに表示させたり、プリントアウトす
ることができることはいうまでもない。
【００６７】デバイス・シミュレーションは半導体デバ
イスの動作についての理解を深めるばかりでなく実際に
半導体デバイスを試作する代りに利用されている。そう
して、デバイス・シミュレーションの結果に基いて、半
導体デバイスの設計を行なうことができる。
イスの動作についての理解を深めるばかりでなく実際に
半導体デバイスを試作する代りに利用されている。そう
して、デバイス・シミュレーションの結果に基いて、半
導体デバイスの設計を行なうことができる。
【００６８】デバイス・シミュレーションを正確かつ短
時間に行えるということは、半導体デバイスの試作費用
の低減や設計を迅速かつ確実に行なうことを可能とする
ということができる。
時間に行えるということは、半導体デバイスの試作費用
の低減や設計を迅速かつ確実に行なうことを可能とする
ということができる。
【００６９】
【発明の効果】以上説明したように本発明は、電界強度
の絶対値がある上限値より大きい高電界領域でのみポテ
ンシャル変動量にダンパーを掛けて反復計算を実行する
ので、基本方程式の非線形性による収束の困難性をなく
すことができかつ不必要に反復回数を増大することを回
避でき計算時間を短縮できるという効果がある。
の絶対値がある上限値より大きい高電界領域でのみポテ
ンシャル変動量にダンパーを掛けて反復計算を実行する
ので、基本方程式の非線形性による収束の困難性をなく
すことができかつ不必要に反復回数を増大することを回
避でき計算時間を短縮できるという効果がある。
【図１】本発明の第１の実施の形態について説明するた
めのフローチャートの部分図である。
めのフローチャートの部分図である。
【図２】図１に続いて示すフローチャートの部分図であ
る。
る。
【図３】移動モデルの一例を示すグラフである。
【図４】本発明の第２の実施の形態について説明するた
めのフローチャートの部分図である。
めのフローチャートの部分図である。
【図５】図４に続いて示すフローチャートの部分図であ
る。
る。
【図６】第１の実施例について説明するためのＭＯＳト
ランジスタの断面図である。
ランジスタの断面図である。
【図７】図６のＭＯＳトランジスタの不純物濃度分布を
示すグラフである。
示すグラフである。
【図８】第１の実施例におけるメッシュを示す図であ
る。
る。
【図９】第１の実施例で使用する電子移動度モデルを示
すグラフである。
すグラフである。
【図１０】第１の実施例で使用する正孔移動度モデルを
示すグラフである。
示すグラフである。
【図１１】第１の実施例におけるポテンシャル変動量の
振舞を示すグラフである。
振舞を示すグラフである。
【図１２】第１の実施例で得られたポテンシャル分布を
示すグラフである。
示すグラフである。
【図１３】第１の実施例で得られた電子移動度の分布を
示すグラフである。
示すグラフである。
【図１４】第２の実施例について説明するためのＭＯＳ
トランジスタの断面図である。
トランジスタの断面図である。
【図１５】図１４のＭＯＳトランジスタの不純物濃度分
布を示すグラフである。
布を示すグラフである。
【図１６】第２の実施例におけるメッシュを示す図であ
る。
る。
【図１７】第２の実施例におけるポテンシャル変動量の
振舞を示すグラフである。
振舞を示すグラフである。
【図１８】第２の実施例で得られたポテンシャル分布を
示すグラフである。
示すグラフである。
【図１９】第２の実施例で得られたインパクトイオン化
率分布を示すグラフである。
率分布を示すグラフである。
【図２０】ニュートン法について説明するためのフロー
チャートである。
チャートである。
【図２１】メッシュの一例を示す図である。
【図２２】カップル法について説明するためのフローチ
ャートである。
ャートである。
【図２３】ガンメル法について説明するためのフローチ
ャートである。
ャートである。
１〜２４ 処理手順 ３１ Ｐ型シリコン基板 ３２ ソース領域 ３３ ドレイン領域 ３４ チャネルドープ領域 ３５ ゲート酸化膜 ３６ 側壁スペーサ １０１〜１０６ 処理手順 １３１〜１３５ 処理手順 １４１〜１４７ 処理手順 Ｂ 基板電極 Ｄ ドレイン電極 Ｇ ゲート電極 Ｓ ソース電極
Claims (3)
- 【請求項１】 半導体デバイスの数値解析のニュートン
反復法による計算の過程で、ｉ回目の計算結果から得ら
れるポテンシャルΨi+1 を用いて電界強度Ｅi+1 を計算
し、前記電界強度の絶対値｜Ｅi+1 ｜が所定の上限値Ｅ
lim を越えた場合、更に電界変動量δＥi ＝Ｅi+1 −Ｅ
i の絶対値｜δＥi ｜が所定の上限値δＥlim を越えた
ときに、ポテンシャル変動量δΨi ＝Ψi+1 −Ψi に所
定のダンパーＤ1iを掛けたＤ1i・δΨi とΨi との和Ψ
i+1 a ＝Ψi ＋Ｄ1i・δΨi をΨi+1 の代りに使用して
（ｉ＋１）回目の計算を行なう手順を有することを特徴
とすデバイス・シミュレーション方法。 - 【請求項２】 半導体デバイスの数値解析のニュートン
反復法による計算の過程で、ｉ回目の計算結果のポテン
シャルΨi を用いて電界強度Ｅi を計算し、前記電界強
度の絶対値｜Ｅi ｜が所定の上限値Ｅlim を越えた場
合、更にｉ回目の計算結果から得られるポテンシャルΨ
i+1 を用いて電界強度Ｅi+1 を計算し、電界強度比の絶
対値｜Ｅi+1 ／Ｅi ｜が所定値Ｒlim を越えたときに、
ポテンシャル変動量δΨi ＝Ψi+1 −Ψi に所定のダン
パーＤ2iをかけたＤ2i・δΨi とΨi との和Ψi+1 b ＝
Ψi ＋Ｄ2i・δΨi をΨi+1 の代りに使用して（ｉ＋
１）回目の計算を行なう手順を有することを特徴とする
デバイス・シミュレーション方法。 - 【請求項３】 電荷保存式、電子電流保存式及び正孔電
流保存式をカップル法で解く請求項１又は２記載のデバ
イス・シミュレーション方法。
Priority Applications (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP7208885A JP2734418B2 (ja) | 1995-08-16 | 1995-08-16 | デバイス・シミュレーション方法 |
| US08/698,375 US5828586A (en) | 1995-08-16 | 1996-08-15 | High speed device simulating method |
| KR1019960033920A KR100219105B1 (ko) | 1995-08-16 | 1996-08-16 | 고속 디바이스 시뮬레이팅 방법 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP7208885A JP2734418B2 (ja) | 1995-08-16 | 1995-08-16 | デバイス・シミュレーション方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPH0955505A JPH0955505A (ja) | 1997-02-25 |
| JP2734418B2 true JP2734418B2 (ja) | 1998-03-30 |
Family
ID=16563738
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP7208885A Expired - Fee Related JP2734418B2 (ja) | 1995-08-16 | 1995-08-16 | デバイス・シミュレーション方法 |
Country Status (3)
| Country | Link |
|---|---|
| US (1) | US5828586A (ja) |
| JP (1) | JP2734418B2 (ja) |
| KR (1) | KR100219105B1 (ja) |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP3885274B2 (ja) * | 1997-03-14 | 2007-02-21 | 富士通株式会社 | シミュレーション用回路変換装置 |
| JP3156764B2 (ja) * | 1997-08-21 | 2001-04-16 | 日本電気株式会社 | 半導体デバイスの衝突電離現象のシミュレーション方法 |
| JP3204300B2 (ja) * | 1997-08-26 | 2001-09-04 | 日本電気株式会社 | デバイスシミュレーション方法 |
| JP3139428B2 (ja) * | 1997-10-31 | 2001-02-26 | 日本電気株式会社 | 拡散シミュレーション方法 |
| JP3382544B2 (ja) | 1998-09-14 | 2003-03-04 | 沖電気工業株式会社 | Ｍｏｓｆｅｔのシミュレーション方法及び装置 |
| JP4281047B2 (ja) * | 2002-03-29 | 2009-06-17 | セイコーエプソン株式会社 | 半導体デバイスシミュレーション方法、半導体素子製造方法 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5684723A (en) * | 1987-11-16 | 1997-11-04 | Fujitsu Limited | Device simulation method and device simulator |
| US5103415A (en) * | 1989-01-13 | 1992-04-07 | Kabushiki Kaisha Toshiba | Computer-simulation technique for numerical analysis of semiconductor devices |
| IE69192B1 (en) * | 1990-12-21 | 1996-08-21 | Hitachi Europ Ltd | A method of generating partial differential equations for simulation a simulation method and a method of generating simulation programs |
| KR100265525B1 (ko) * | 1994-09-09 | 2000-09-15 | 가네꼬 히사시 | 반도체 디바이스 시뮬레이션의 초기 포텐셜 값 추정 방법 |
-
1995
- 1995-08-16 JP JP7208885A patent/JP2734418B2/ja not_active Expired - Fee Related
-
1996
- 1996-08-15 US US08/698,375 patent/US5828586A/en not_active Expired - Fee Related
- 1996-08-16 KR KR1019960033920A patent/KR100219105B1/ko not_active IP Right Cessation
Also Published As
| Publication number | Publication date |
|---|---|
| KR970013164A (ko) | 1997-03-29 |
| US5828586A (en) | 1998-10-27 |
| JPH0955505A (ja) | 1997-02-25 |
| KR100219105B1 (ko) | 1999-10-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Chevalier et al. | Si/SiGe: C and InP/GaAsSb heterojunction bipolar transistors for THz applications | |
| JP2734418B2 (ja) | デバイス・シミュレーション方法 | |
| CN106339562B (zh) | 具有电流守恒特性的能带间非局域量子隧穿模拟方法 | |
| Carroll et al. | Degradation of silicon bipolar junction transistors at high forward current densities | |
| Gildenblat et al. | Introduction to PSP MOSFET model | |
| Ghosh et al. | Analysis of kink reduction and reliability issues in low‐voltage DTD‐based SOI TFET | |
| Hareland et al. | A simple model for quantum mechanical effects in hole inversion layers in silicon PMOS devices | |
| Lee et al. | Numerical modeling of advanced semiconductor devices | |
| Ghetti et al. | Self-consistent simulation of quantization effects and tunneling current in ultra-thin gate oxide MOS devices | |
| JUNGEMANN et al. | Full-band Monte Carlo device simulation of a Si/SiGe-HBT with a realistic Ge profile | |
| Zekry et al. | A comprehensive semi-analytical model of the polysilicon emitter contact in bipolar transistors | |
| US20020010564A1 (en) | Semiconductor device simulation method | |
| US11742247B2 (en) | Epitaxial growth of source and drain materials in a complementary field effect transistor (CFET) | |
| Daliento et al. | An improved model for extraction of strongly spatial dependent lifetimes with the ac lifetime profiling technique | |
| Yu et al. | Equivalent circuit approach for single electron transistor model for efficient circuit simulation by SPICE | |
| Bose et al. | Analytical model and simulation‐based analysis of a work function engineered triple metal tunnel field‐effect transistor device showing excellent device performance | |
| JP3221354B2 (ja) | デバイスシミュレーション方法 | |
| Mattausch et al. | HiSIM2 Circuit simulation-Solving the speed versus accuracy crisis | |
| JP3926150B2 (ja) | デバイス・シミュレーション方法およびデバイス・シミュレーション装置 | |
| Palestri et al. | A drift-diffusion/Monte Carlo simulation methodology for Si/sub 1-x/Ge/sub x/HBT design | |
| JP3204300B2 (ja) | デバイスシミュレーション方法 | |
| Selvi et al. | Potential model and simulation analysis of dual material gate vertical TFET with impact of interface trap charges | |
| Ferguson et al. | An approximate two-dimensional transient analysis of charge transport in a buried-channel charge-coupled device | |
| JP3214108B2 (ja) | 半導体デバイスのシミュレーション方法 | |
| JP3152184B2 (ja) | デバイスシミュレーション方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 19971125 |
|
| LAPS | Cancellation because of no payment of annual fees |