JPH03253056A - 半導体デバイスの動作解析法、所定の物理現象の解析法およびそのために使用する装置 - Google Patents

Abstract

(57)【要約】本公報は電子出願前の出願データであるた
め要約のデータは記録されません。

Description

【発明の詳細な説明】 〔発明の目的〕 （産業上の利用分野） 本発明は、半導体デバイス及び所定物理現象のモデリン
グに関する。

（従来の技術） 半導体デバイスモデリングにおいて用いられる基本方程
式は、通常法の形で表わされる。以下の式において、τ
は物理時間を表わす。

ａｐ／δτ＝（１／ｑ）　（ｄｉｖ　Ｊｐ）＋Ｇ　　Ｕ
　　・＝　Ｑ）ａｎ／∂ｔ＝（１／ｑ）（ｄｉｖＪｎ）
十〇−Ｕ　　　＝−■Ｊｐ＝−ｑｄｐ（ｇｒａｄｐ）−
ｑμｐｐ（ｇｒａｄψ）・・・■Ｊｎ＝−ｑｄｎ（ｇｒ
ａｄｎ）−ｑμｎｎ（ｇｒａｄψ）　　　−Ｋｄｉｖ　
（ｇｒａｄψ）＝（−ｑ／１）（Ｎｄ−Ｎａ十ｐ−ｎ）
−０５１これらの方程式は、従来法においても、本発明
をその一部として含む方法においても共通に用いられる
ものである０式（１）および■は各々正孔、電子の連続
方程式を表わす。Ｇ、Ｕは各々キャリアの発生、消滅を
表わす。式（３）、■は各々正孔電流密度、電子電流密
度を具体的に記述したものであり、正孔電流密度ＪＰお
よび電子電流密度Ｊ。が、それぞれキャリア密度の傾斜
ｇｒａｄ　ｐ　、　ｇｒａｄ　ｎに比例する拡散電流成
分（第１項）と、電位傾斜ｇｒａｄφと各キャリア密度
に比例するドリフト成分（第２項）の代数和からなるこ
とを表わす。弐〇は、電位ψと空間電荷密度の関係を表
わすポアソン方程式である。

以下、従来の解法について具体的に説明する。

方程式（１）〜■の数値解を求めるための第１のステッ
プは、これらの式を連続形（微分方程式）から離散系（
差分方程式）に書き換えることである。

この際、解析の対象となるデバイス空間を長方形のメツ
シュに分割する必要がある。二次元空間モデルについて
メツシュ分割を概念的に第１図に示す。縦、横の分割点
間隔ｈｘ、ｈｙを第１図のように定義し、かつ、各々の
中間点の距離ｈＸ、ｈｙを次式で定義する。

ｈＸ’（Ｍ）＝（１／２）［ｈＸ（Ｍ’−１）＋ｈＸ（
Ｍ’）］−（］６−１ｈｙ’（Ｎ）＝（１／２）［ｈｙ
（Ｎ’　−１）＋ｈｙ（Ｎ’）コ・・　　（６−２）式
（３）、　（６−１）、　（６−２）より、離散系では
、正孔電流の発散ｄｉｖＪｐは次の近似式で表わせる。

ｄｉｖ　Ｊｐ”ａ　Ｊｐｘ／　ａ　ｘ＋ａ　Ｊｐｘ／　
ａ　Ｙ＝　　［ＪＰＸ（Ｍ’）−ＪＰＸ（Ｍ’−１）コ
／ｈＸ′（Ｍ）＋［Ｊｐｙ（Ｎ’）　　Ｊｐｙ（Ｎ’　
　１）］／ｈｙ’（Ｎ）・・・ω 電子電流についても同様の書き替えが可能であるが、自
明であるので省略する。

以上の式により、式（１）、■全体を書き替えた結果は
次の通りである。

（１／（１）（［ＪＰＸ（Ｍ’）−ＪＰＸ（Ｍ’−１）
／ｈｘ’（Ｍ）＋［Ｊｐｙ（Ｎ’）−Ｊｐｙ（Ｎ’−１
）］／ｈｙ’（Ｎ））＝−ｕ（Ｍ、Ｎ）・・・（８−１
） （１／　ｑ）（［Ｊｎｘ（Ｍ’）−Ｊｎｘ（Ｍ’　　１
）／　ｈｘ’（Ｍ）＋［Ｊｎｙ（Ｎ’）−Ｊｒ＋ｙ（Ｎ
’−ｉ）コ／　ｈｙ’　（Ｎ））＝　　　Ｕ（Ｍ、Ｎ）
・・・（８−２） ただし、（８−１）、（８−２）式において、キャリア
発生項Ｇは通常のデバイス解析では必要とされないので
、簡単化のため省略した。また、同じく簡単化のため、
以下においては直流定常解の計算を行う場合を考えるも
のとし、物理的時間微分項ａｐ／ａτおよびａ　ｎ　／
　ａτはゼロとおいた。

しかし、これらの物理的時間微分項がゼロでない非定常
問題の場合についても、本質的な差異を生じることなく
計算を実行できることを指摘しておく。

同様にしてポアソン方程式■は、次のように表わされる
。

［ｔ／ｈｘ’（Ｍ）］　（［ψ（Ｎ＋１．Ｎ）−ψ（Ｍ
、　Ｎ）Ｍｈｘ（Ｍ’）−ψ（Ｍ、　Ｎ）−ψ（Ｍ−１
、Ｎ）］／　ｈ　Ｘ（Ｍ’−１））＋［１／ｈｙ’（Ｎ
）］　（［ψ（Ｍ、Ｎ＋１）−ψ（Ｍ、　Ｎ）］／ｈｙ
（Ｎ’）−［ψ（Ｍ、Ｎ）−ψ（Ｍ、　Ｎ−１）］／ｈ
ｙ（Ｎ’−ｘ））＝　　　（ｑ／ε）［ｒ（Ｍ、　　Ｎ
）＋ｐ（Ｍ、　　Ｎ）−ｎ（Ｍ、　　Ｎ）コ・・・■ ただし、簡単のため、ドナー、アクセプタ不純物濃度の
差をｒ＝Ｎｄ−Ｎａとおいた。

ここで、以下の展開のために、キャリア再結合項Ｕ　（
Ｍ、Ｎ）を具体的にＳ　ｈｏｃｋｌｅｙ−Ｒｅａｄ−Ｈ
ａｌｌの形で与える。即ち。

Ｕ（Ｍ、　Ｎ）＝ ［ｐ　（Ｍ、Ｎ）・ｎ　（Ｍ、Ｎ）−ｎ　ｔ”（Ｍ、Ｎ
）］／　（τｎ［ｐ　（Ｍ、Ｎ）＋ｎ４（Ｍ、Ｎ）コ＋
でｐ［ｎ（Ｍ、Ｎ）＋ｎ４（Ｍ、Ｎ）］・・・　（１０
） 次に、電流方程式■、■を離散化する。この際、従来か
ら知られているように、数値計算上の安定性を確保する
ために、Ｓ　ｃｈａｒｆｅｔｔｅｒ−Ｇ　ｕｗａｍｅｌ
の近似形を採用する。その結果、正孔、電子の電流密度
のｘ、ｙ成分は各々次のように与えられる。

Ｊｐｘ（Ｍ’）＝［（１／　ｈｘ（Ｍ’）］Ｘ［λＰＸ
１（Ｍ’）Ｐ　（ＭｔＮ）＋λｐｘ２（Ｍ’）Ｐ　（Ｍ
＋　１　、Ｎ）］　　　　　　　・＝　（１１−１）Ｊ
　ｐｙ（Ｎ’　）”　［ｑ／　ｈ　ｙ（Ｎ’　）］　Ｘ
　［λｐｙ□（Ｎ’）ｐ　（Ｍ、Ｎ）＋λｐｙｚ（Ｎ’
）ｐ（Ｍ＋Ｎ＋１）］　　　　　　　・・・（１１−２
）Ｊｎｘ（Ｍ’）＝［ｑ／　ｈｘ（Ｍ’）］Ｘ［λｎｘ
ｘ（Ｍ’）ｎ　（Ｍ、Ｎ）＋λｎｘ□（Ｍ’）ｎ　（Ｍ
＋　１　、Ｎ）］　　　　　　　−（１１−３）Ｊｎｙ
（Ｎ’）＝［ｑ／ｈｙ（Ｎ’）］Ｘ［λｎｙｌ　（Ｎ’
）　ｎ　（ＭＩＮ）十λｎｙ２（Ｎ’）ｎ（Ｍ十Ｎ＋１
）］　　　　　　　　”’　（１１−４）ただし、上式
のλは次式で与えられる。

λｐｘ、（Ｍ’）＝［μＰ（Ｍ’）／θ（Ｍ’）］Ｘ［
β（Ｍ’）／（１−ｅ−β（ｆ））］λｐｘｚ　（Ｍ’
　）＝［μＰ（Ｍ’）／θ（Ｍ’）］Ｘ［β（Ｍ’）／
（１−ｅ−”’−’）］λｐｙｔ（Ｎ’）＝［／’ｐ（
Ｎ’）／θ（Ｎ’）］Ｘ［β（Ｎ’）／（１−ｅ″″β
（Ｎ−））］λＰｙ２（Ｎ′）＝［μＰ（Ｎ′）／θ（
Ｎ′）コ×［β（Ｎ　′）／（ｌ、−ｅａ（Ｎ　　））
］λｎＸ、（Ｍ’）＝［Ｕｎ（Ｍ′）／μＰ（Ｍ’）］
Ｘ［λＰＸ２　（Ｍ’　）λｎｘ２（Ｍ’）＝［Ｕｎ（
Ｍ′）／μＰ（Ｍ’）コ×［λｐｘ−（Ｍ’）λｎｙ、
（Ｎ’）＝［Ｕｎ（Ｎ′）／μｐ（Ｎ’）コ×［λｐｙ
−（Ｎ’）λｎｙ、（Ｎ’）＝［Ｕｎ（Ｎ′）／μｐ（
Ｎ’）］Ｘ［λｐｙ１（Ｎ’）・・・（１２） 以上により、原方程式（１）〜■の離散系への書き替え
が終了した。即ち、我々が解くべき方程式は、（８−１
）、（８−２）式に電流の式（１１−１）ないしく１ｌ
−４）、補助関係式（１２）と再結合の式（１０）とを
代入したものと、式■）であり、これらは３個の未知量
ｐ、ｎ、ψを有する３個の方程式からなる。

ところが、（８−１）、（８−２）式を見ると、これら
は未知量につき非線形な関係式を含むから、このままで
は解を求めることが困難である。そこで、次に、これら
の式を未知量につきティラー展開し、２次以上の項を無
視することにより、未知量の変化分δｐ、δｎ、δψに
ついての線形方程式を導出する。

まず、各未知量を次のように、試行値（上付きのゼロで
表わす）と、変化分の和からなるものとする。

ｐ（Ｍ＋　Ｎ）＝ｐ’（Ｍｙ　Ｎ）＋δｐ（Ｍ、Ｎ）ｎ
（Ｍ、Ｎ）＝ｎ’（Ｍ、Ｎ）＋δｎ（Ｍ、Ｎ）ψ（Ｍ、
Ｎ）＝ψ’（Ｍ、Ｎ）十δψ（Ｍ、Ｎ）・・・（１３） そこで電流の式（１１−１）　〜（１１−４）をＴａｙ
ｌｏｒ展開し、２次以上の項を無視すると、次の結果を
得る。

ＪＰＸ（Ｍ’）〜、Ｔｐｘ。（Ｍ′） ＋［ａＪｐｘ。（Ｍ’）／ａｐ（Ｍ、Ｎ）］Ｘδｐ（Ｍ
、Ｎ）＋［ａＪｐｘ。（Ｍ’）／　ａ　ｐ　（Ｎ＋１　
、Ｎ）コ×δ　ｐ（Ｎ＋１．Ｎ）十［ａ　Ｊ　ｐｘ”　
（Ｍ’）／　ａ　ψ（Ｍ、Ｎ）ｌｘ　８　ψ（Ｍ、Ｎ）
＋［ａＪｐｘ。（Ｍ’）／ａ　ψ（Ｎ＋１　、Ｎ）］Ｘ
δψ（Ｎ＋１．Ｎ）・・・（１４−１） ＪＰｙ（Ｎ’）〜ＪＰｙ。（Ｎ′） ＋［ａＪｐｙ。（Ｎ’）／　ａ　ｐ　（Ｍ、Ｎ）］Ｘδ
ｐ（Ｍ、Ｎ）十［ａＪｐｙ。（Ｎ’）／　ａ　ｐ　（Ｍ
、Ｎ＋１　）］Ｘ　８　ｐ　（Ｍ、Ｎ＋１）＋［ａ　Ｊ
ｐｙ。ＣＮ’）／ａ　　ψ（Ｍ、Ｎ）コＸ　８　　ψ（
Ｍ、Ｎ）＋［ａＪｐｙ。（Ｎ’）／ａ　ψ（Ｍ、Ｎ＋１
）］Ｘδψ（Ｍ、Ｎ＋１）・・・（１４−２） ＪｎＸ（Ｍ′）〜ＪｎＸ。（Ｍ′） ＋［ａＪｎｘ。（Ｍ’）／ａｎ（Ｍ、Ｎ）］ｘδｎ（Ｍ
、Ｎ）＋　［ａ　Ｊ　ｎｘ。（Ｍ’　）／　ａ　ｎ　（
Ｍ＋　１　、　Ｎ　）］　ｘδｎ（Ｎ＋１．Ｎ）＋［ａ
Ｊｏｘ。（Ｍ’）／ａψ（Ｍ、Ｎ）］Ｘδψ（Ｍ、Ｎ）
＋［ａＪｎｘ。（Ｍ’）／ａ　ψ（Ｎ＋１　、Ｎ）コ×
δ　ψ（Ｎ＋１．Ｎ）・・・（１４−３） Ｊｎｙ（Ｎ′）〜Ｊｎｙ。（Ｎ′） ＋［ａＪｎｙ。（Ｎ’）／　ａ　ｎ　（Ｍ、Ｎ）］Ｘδ
ｎ（Ｍ、Ｎ）＋［ａ　Ｊ　ｎｙ。（Ｎ’）／ａｎ（Ｍ、
Ｎ＋１）］ｘδ　ｎ（Ｍ、Ｎ＋１）＋［ａＪｎｙ。（Ｎ
’）／　ａ　　ψ（Ｍ、Ｎ）コ×δ　ψ（Ｍ、Ｎ）＋［
ａＪｎｙ。（Ｎ’）／δψ（Ｍ、Ｎ＋１）ＩＸδψ（Ｍ
、Ｎ＋１）・・・（１４−４） ただし上つきゼロをもっ容量は未知量の試行値Ｐ’ｔ　
ｎ”、ψ０に対応するものとする。

電流のほか再結合項Ｕ　（Ｍ、Ｎ）も非線形なので、こ
れも同様の計算を行うと、つぎの結果を得る。

Ｕ（Ｍ、Ｎ）＝Ｕ’（Ｍ、Ｎ）＋ＵＰ。（Ｍ、Ｎ）Ｘδ
ｐ　（Ｍ　ｔ　Ｎ　）　＋　ｕ　ｎ。（Ｍ、Ｎ）δｎ（
Ｍ、Ｎ）ただし Ｕｐ（Ｍ、Ｎ）＝ａ　Ｕ（Ｍ、Ｎ）／δｐ（Ｍ、Ｎ）＝
［ｎ　（Ｍ、Ｎ）　　ｘ　ｎＵ（Ｍ、Ｎ）／　（ｘ。

Ｘ［ｐ　（Ｍ、Ｎ）＋　ｎｉ（Ｍ、Ｎ）］＋　τｐ［ｎ
　（Ｍ、Ｎ）＋ｎ　ｉ（Ｍ　−Ｎ　）］ Ｕｎ（Ｍ　＝　Ｎ　）　＝ａ　Ｕ　（Ｍ　−Ｎ）／　ａ
　ｎ　（Ｍ　−Ｎ）＝［ｐ（Ｍ、Ｎ）−τ、Ｕ（Ｍ、Ｎ
）／（τ。

Ｘ［ｐ　（Ｍ、Ｎ）＋ｎｉ（Ｍ、Ｎ）コ＋τｐ［ｎ（Ｍ
、Ｎ）＋ｎｌ（Ｍ、Ｎ）コ）　　　　　　　　　　　　
　　　　　　　　−（ｔＳ）ここで、（１３）、（１４
−１）、（１４−２）、（１４−３）。

（１４−４）、　（１５）を（８−１）、（８−２）お
よび■）に代入し、項を整理すると、っぎの行列・ベク
トル方程式が得られる。

ＡＴ　（Ｍ、Ｎ）ｅＴ（Ｍ、Ｎ−１） 十ＢＴ　（Ｍ、　Ｎ）　ｅｔ　（Ｍ、　Ｎ）＋ＣＴ　（
Ｍ、Ｎ）　　ｅＴ　（Ｍ、Ｎ＋１）＋ＤＴ　（Ｍ、Ｎ）
ｅｒ　（Ｍ　　１．Ｎ）＋ＥＴ　（Ｍ、Ｎ）ＯＴ（Ｎ＋
１．Ｎ）＝ＦＴ（Ｍ、Ｎ） ただし ・・・　（１６） ・・・　（１７） なお、各行列、ベクトル要素を具体的に書き下すと、つ
ぎのようになる。

ａ□ｘ＝　　［１／　ｑ　ｈｙ’（Ｎ）］ＸａＪｐｙ。

（Ｎ−１）／ａｐ（Ｍ、　Ｎ−１）ａ１□＝Ｏ ａｔａ＝−［１／ｑ　ｈｙ’（Ｎ）コ ×δＪ　ｐｙ’　（Ｎ　　１　）／δψ（Ｍ、Ｎ−１）
ａ２１＝Ｏ ａ２□＝　　［１／　ｑ　ｈｙ’（Ｎ）］ＸａＪｎｙ。

（Ｎ−１）／ａｎ（Ｍ、　Ｎ−１）ａｚ３＝　　［１／
ｑｈｙ’（Ｎ）］ ×δＪｎｙ。（Ｎ−１）／ａψ（Ｍ、Ｎ−１）ａ、、＝
Ｏ ａ３２＝Ｏ ａａ３＝”ＩｘＣＭ＋　　Ｎ） ＝１／［ｈｙ’（Ｎ）ｂｙ（Ｎ−１）］ｂ１□＝［１／
ｈＸ′（Ｍ）］ ×［δＪ　ＰＸ”　（Ｍ）／　ａ　ｐ　（Ｍ、　Ｎ）−
δ、Ｊｐｘ。（Ｍ−１）／　ａ　Ｐ　（Ｍ、　Ｎ）］＋
　［１／　ｈ　ｙ’　（Ｍ）］ ｘ［ａ　Ｊｐｙ’ＣＭ）／δｐ（Ｍ、Ｎ）−ａ　Ｊ　ｐ
ｙ’　（Ｍ−１）／　ａ　ｐ　（Ｍ、　Ｎ）］−１−ａ
Ｕ’（Ｍ、Ｎ）／ａｐ（Ｍ、Ｎ）ｂ、、＝ａＵ’ＣＭ、
Ｎ）／ａｎ（Ｍ、Ｎ）ｂ　１！　＝　［１／　ｑｈ　ｘ
’　（Ｍ）］ｘ［ａＪｐｘ。（Ｍ）／ａ　ｌ／ｌ　（Ｍ
、　Ｎ）−ａＪＰＸ。（Ｍ−１）／　ａ　　ψ（Ｍ、　
　Ｎ）］＋［１／　ｑ　ｈｙ’（Ｍ）コ Ｘ［ａＪｐｙ。（Ｍ）／ａ　ψ（Ｍ、Ｎ）ａＪｐｙ。Ｃ
Ｍ−１）／ａ　ψ（Ｍ、Ｎ）］ｂ２１＝−ａＵ’（Ｍ、
Ｎ）／ａｐ（Ｍ、Ｎ）ｂ２□＝　［１／　ｑ　ｈ　ｘ’
　（Ｍ）］Ｘ　［ａ　Ｊ　ｎｘ’　（Ｍ）／　ａ　ｎ　
（Ｍ、　Ｎ）ａ　Ｊ　ｎｘ’　（Ｍ　　１　）／　ａ　
ｎ　（Ｍ、　Ｎ）］＋［１／　ｑ　ｈ　ｙ’　（Ｍ）］ Ｘ　［ａ　Ｊ　ｎｙ’　ＣＭ）／　ａ　ｎ　（Ｍ、　Ｎ
）ａ　Ｊ　ｎｙ’　（Ｍ　　　１　）／　ａ　ｎ　（Ｍ
、　　Ｎ）コ−ａＵ’（Ｍ、Ｎ）／ａｎ（Ｍ、Ｎ） ｂ　２３　＝　［１／　ｑｈ　ｘ’　（Ｍ）］Ｘ［ａＪ
ｎｙ’（Ｍ）／ａ　ψ（Ｍ、Ｎ）ａ　Ｊ　ｎｘ’　（Ｍ
　　　ｌ　）／　ａ　４’　（Ｍ、　Ｎ）コ十［１／　
ｑ　ｈ　ｙ’　（Ｍ）］ ｘ［ａＪｎｙ’ＣＭ）／ａ　ψ（Ｍ、Ｎ）ａ　Ｊ　ｎｙ
’　（Ｍ　　１　）／δψ（Ｍ、Ｎ）］ｂ３１＝ｑ／ε
　　　　ｔ　　ｂｉ□＝−ｑ／εｂ３３＝　［１／ｈｘ
’（Ｍ）］［１／ｈｘ（Ｍ）＋１／ｈｘ（Ｍ−ｉ）コ −［１／ｈｙ’（Ｎ）＋１／ｈｙ（Ｎ−１）］Ｃ□１＝
［１／　ｑ　ｈｙ’（Ｎ）コ ×［δＪＰ３゜（Ｎ）／　ａ　ｐ　（Ｍ、　Ｎ　＋　１
　）］Ｃ□よ＝Ｏ ｃ　ｚ３　＝　［１／　ｑ　ｈ　ｙ’　（Ｎ）］Ｘ［ａ
Ｊｐｙ。（Ｎ）／ａ　ψ（Ｍ、Ｎ＋　１）］ｃ２□＝Ｏ ｃ２□＝［↓／　ｑ　ｈ　ｙ’　（Ｎ）］Ｘ［ａＪｎｙ
。（Ｎ）／　ａ　　ｎ　（Ｍ、　　Ｎ＋　１）コｃ　ｚ
３＝　［１／　ｑ　ｈ　ｙ’　（Ｎ）］Ｘ［ａ　Ｊ　ｎ
ｙ’（Ｎ）／　ａ　φ（Ｍ、　Ｎ＋　１）］Ｃ３ｘ　＝
　Ｏ＋　　Ｃａ　ｚ　＝　Ｏｏ　３ｍ　＝　１　／　ｈ
　ｙ’　（Ｎ）　ｈ　ｙ（Ｎ）ｄ１□＝　　［１／ｑｈ
ｘ’（Ｍ）］ Ｘ［ａＪｐｘ。（Ｍ−１）／　ａ　　ｐ　（Ｍ−１、Ｎ
）コｄｉ２＝Ｏ ｄ　１ｚ　＝　　［１／　ｑｈ　ｘ’　（Ｍ）］×［ａ
ＪＰＸ。（Ｍ−１）／ｃｌ　ψ（Ｍ−ｌ、Ｎ）］ｄ２□
二〇 ｄ２□＝　　［１／ｑｈｘ’（Ｍ）］ Ｘ［ａＪｎｘ’（Ｍ　　　１）／ａｎ（Ｍ　　　１　、
　Ｎ）コｄｚｉ＝−［１／ｑｈｘ’（Ｍ）］ ＸＣａＪｎｘ’（Ｍ　　１）／δψ（Ｍ−１，Ｎ）］ｄ
３．＝Ｏ、ｄ、ｚ＝０ ｄ３ｉ＝１／ｈｘ’（Ｍ）ｈｘ（Ｍ−１，）ｅ　、ｘ　
＝　［１／　ｑ　ｈ　ｘ’　（Ｍ）］Ｘ［ａ　Ｊｐｘ。

ＣＭ）／　ａ　ｐ　（Ｍ＋　１　、　　Ｎ）］ｅ１□＝
Ｏ ｅ□３　＝　［１／　ｑ　ｈ　ｘ’　（Ｍ）コＸ［ａ　
Ｊｐｘ。（Ｍ）／　ａ　　ψ（Ｍ＋　１　、　　Ｎ）］
ｅ２１＝：Ｏ ｅ、２＝［１／　ｑ　ｈｘ’（Ｍ）コ Ｘ［ａＪｎｘ。（Ｍ）／　ａ　ｎ　（Ｍ＋　１　、　Ｎ
）］ｅ２１Ｉ＝［１／ｑｈｘ’（Ｍ）］ ×［ａＪ　ｏｘ’　（Ｍ）／δψ（Ｍ＋１．Ｎ）１６３
１”Ｏｔ　　Ｂ３２”Ｏ ｅ　３３＝　１　／　ｈ　ｘ’　（Ｍ）　ｈ　ｘ（Ｍ）
ｆ工＝　　Ｕ’（Ｍ、Ｎ）　十Ｕｐ。（Ｍ、Ｎ）ｐ’Ｃ
Ｍ、Ｎ）十Ｕ♂（Ｍ、Ｎ）ｎ’（Ｍ、Ｎ） ［１／　ｑ　ｈ　ｘ’　（Ｍ）］ ×［ａＪＰ。（Ｍ−ＬＶａ　ψ（Ｍ−１，Ｎ）］ＸＩＪ
’（Ｍ−１，Ｎ）−φｏ（Ｍ、Ｎ）］十［１／　ｑ　ｈ
ｘ’（Ｍ）コ Ｘ［ａ　ＪＰ。（Ｍ）／ａ　ψ（Ｍ、Ｎ）コ×［ψ’（
Ｍ、　　Ｎ）−ψ’（Ｍ＋１．Ｎ）コ［１／ｑｈｙ’（
Ｎ）］ Ｘ［ａＪ、。（Ｎ−１）／ａ（＃（Ｍ、Ｎ−１）］×［
ψ’（Ｍ、　　Ｎ−１）−ψ’（Ｍ、　　Ｎ）コ＋［１
／　ｑ　ｈｙ’（Ｎ）コ ×［ａＪＰ。（Ｎ）／ａψ（Ｍ、Ｎ）］Ｘ［ψ’（Ｍ、
　　Ｎ）　−ψ’（Ｍ、　　Ｎ＋１）］ｆ２＝Ｕａ（Ｍ
、　Ｎ）−ＵＰ。（Ｍ、　Ｎ）ｐ’（Ｍ、　Ｎ）−Ｕｎ
。（Ｍ、Ｎ）ｎ’（Ｍ、Ｎ） ［１／　ｑ　ｈ　ｘ’　（Ｍ）］ Ｘ［ａ　Ｊ、’（Ｎ−１）／δ　ψ（Ｍ−１，Ｎ）コ×
［φ’（Ｍ−１，Ｎ）−ψ’（Ｍ、Ｎ）］十［１／　ｑ
　ｈ　ｘ’　（Ｍ　）］ ×［ａＪｎ。（Ｍ）／ａ　ψ（Ｍ、　　Ｎ）コｘ［φ’
ＣＭ、　　Ｎ）　−ψＱ（Ｍ＋　１．　　Ｎ）］［１／
　ｑ　ｈｙ’（Ｎ）コ ×［ａＪ♂（Ｎ−１）／ａψ（Ｍ、Ｎ−１）］Ｘ［ψ’
（Ｍ、Ｎ−１）−ψ’（Ｍ、Ｎ）］十［１／　ｑ　ｈ　
ｙ’　（Ｎ）］ ｘ　［ａＪ　ｎ″（Ｎ）／　ａ　ψ（Ｍ、　Ｎ）］×［
ψ’（Ｍ、　　Ｎ）−ψｆｌ（Ｍ、　　Ｎ＋１）コＬ＝
　　（ｑ／ε）Ｉ”（Ｍ、Ｎ） ・・・（１８） 以上により、離散化した非線形方程式（８−１）。

（８−２）、■を線形化した結果、行列・ベクトル方程
式（１６）が導かれることがわかった。

式（１６）は、２次元配列された多数の分割点の１個分
に相当するものであるから、デバイス平面全体について
解を得るためには式（１６）を全分割点につき連立した
ものを解かねばならない。いま仮にデバイス全体に対応
する分割点番号が、工≦Ｍ≦に、１≦Ｎ≦Ｌの範囲にあ
るものとすれば、式（１６）を全体につき書き下した結
果は第２図に示すようになる。但しここで簡単のため、
Ａ丁（Ｍ、Ｎ）をＡＮＮ、　ｅ（Ｍ　、　Ｎ　）をｅＭ
Ｎ＋　Ｅ丁ＣＭＩ　Ｎ）をＦＭＮなどに書きかえた。

第２図において、Ｂ１１など個々の行列は（３×３）の
ディメンションをもつから、全体の行列のディメンショ
ンは（３ＫＬ）ｘ　（３ＫＬ）となる。

いま標準的なデバイスモデルとしてに＝３０．Ｌ＝４０
を与えるものとすれば、全体の行列のディメンションは
３，６００　Ｘ　３，６００　＝　１２，９６０，００
０　＝　１．２９６　Ｘ　１０７の大きさとなる。

第２図の行列・ベクトル方程式を実際に解く方法は色々
あり、最も普通に用いられる方法はＧ　ａｕｓｓの消去
法である。今日計算機の発達により大容量メモリの高速
計算が可能となったため、上記の数値例の程度の大きさ
の問題は消去法を用いて解くことができる。別な方法と
しては各種の反復法があり、これは行列サイズが非常に
大きい場合、有効な手段となる。

行列解法はデバイスモデリングに限らず、他の物理的、
光学的問題に広く応用されており、既によく知られた計
算プロセスであるので、以下その解は所定の計算時間を
もって求まったものとする。

すなわち式（１６）の解が、与えられたデバイスの空間
全体につき求まったものとする。このことは、非線形方
程式（８−１）、　（８−２）および■を線形近似した
ものの解が求まったことに他ならないゎけであるが、こ
の解を非線形方程式に代入してみると、非線形効果があ
るため、一般にその両辺は等しくない。そこで既に求ま
った［ｐ（Ｍ、Ｎ）。

ｎ（Ｍ　ｇ　Ｎ　）　ｔ　ψ（Ｍ、Ｎ）、１≦Ｍ≦に、
１≦Ｎ≦Ｌコを改めて［：ｐ’（Ｍ＋　Ｎ）ｒ　ｎ’（
Ｍ、Ｎ）ｖψ’（Ｍ、Ｎ）、１≦Ｍ≦に、１≦Ｎ≦Ｌコ
と置き換え、再度上述したのと同一の計算プロセスを実
行する。以下同様にして、必要回数だけ反復計算を実行
すれば、遂には非線形方程式（８−１）、（８−２）、
■を完全に満足する解が求まる。

（発明が解決しようとする課題） 上述した行列解法は、今日の計算機の性能向上により、
メモリの大容量化と高速計算が可能となったため、かな
り大規模な問題に対してもこれを適用することができる
。しかし計算機の性能向上に伴って対象となる問題がさ
らに大形化するため、行列解法の高速化はモデリングの
実行上、依然として最大問題のひとつである。

この問題の解法のため１行列解法を適用することなしに
所定の方程式の解を求め、半導体デバイスの動作解析、
所定の物理現象の解析を行う方法及び装置が提案されて
いるが、本発明は、そこにおける空間電荷の関係又は増
幅利得又は係数行列の決定法を改良したものである。

〔発明の構成〕

（課題を解決するための手段） 本発明の第１の半導体デバイスの動作解析法における空
間電荷の関係決定法は、半導体デバイスモデリングにお
ける。電位と空間電荷の関係を表わすポアソン方程式を
コンピュータにより解析する方法において、前記ポアソ
ン方程式を時間微分項ａφ／∂ｔと空間電荷の関係λを
含む下記（ａ２）式に置き換え、 ∂ψ／∂ｔ＝λｆ＋　　　　　　　　−（ａ２）半導体
デバイスの分割直配ｍｌ　（Ｍ、　Ｎ）を決定して上記
（ａ２）式を下記の（ｂ２）式に変換し、ａ　ψＣＭ、
　Ｎ）／ａ　ｔ＝λ（Ｍ、　Ｎ）ｆ＋（Ｍ、　Ｎ）・・
・（ｂ２） 上記（ｂ２）式を時間積分することにより、前記ポアソ
ン方程式の解を求める動作解析法において、空間電荷の
関係λを、安定性条件から決定することを特徴とする。

本発明の第２の半導体デバイスの動作解析法における空
間電荷の関係決定法は、半導体デバイスモデリングにお
ける、電子、正孔の輸送方程式およびポアソン方程式か
らなる連立方程式系をコンピュータを用いて解析する方
法において、前記連立方程式系　ｆ　ｐ：　Ｏ、ｆ　ｎ
”　Ｏｔ　ｆ　＊＝　Ｏを、人口的時間微分項ｄ　ｐ　
／　ｄ　ｔ　／　ｓ　ｄ　ｎ　／　ｄ　ｔ　＝　ｄψ／
ｄｔと空間電荷の関係λ、λ。、λφとを含む下記（ａ
　１）、　（ｂｌ）、（ｃ　１）式に置き換え、ｄ　ｐ
／ｄ　ｔ＝−λＰｆｐ・・・（ａｌ）ｄｎ／ｄｔ＝　　
λｎｆｎ　　　　　　　・（ｂ　ｌ）ｄψ／ｄｔ＝　　
λ＄ｆ＊　　　　　　　・＝　（ｃ　１）半導体デバイ
スの分割点配置（Ｍ、Ｎ）を決定し、上記（ａ　１）、
　（ｂ　１）、　（ｃ　１）式を下記の（ｄ）。

（ｅ）、（ｆ）式に変換し、 ｄ　ｐ　（Ｍ、Ｎ）／　ｄ　ｔ　＝−λｐ　（Ｍ　、　
Ｎ　）　ｆ　ｐ　（Ｍ　、　Ｎ　）・・・（ｄ） ｄ　ｎ　（Ｍ、Ｎ）／　ｄ　ｔ　＝λ、（Ｍ、Ｎ）ｆ、
（Ｍ、Ｎ）（ｅ） ｄψ（Ｍ、Ｎ）／ｄ　ｔ　＝λ◆（Ｍ、Ｎ）ｆ◆（Ｍ、
Ｎ）・・・　（ｆ） は次式のと ただし各式の右辺のｆＰ＋　ｆｎ＋　ｆ◆おり定義し、 ｆ　ｐ（Ｍ、Ｎ）＝（１／　ｑ）　（［Ｊｐｘ（Ｍ）Ｊ
　ｐｘ（Ｍ　　１　）］／　ｈ　ｘ’　（Ｍ））＋　（
１／　ｑ）　（［Ｊ　ｐｙ（Ｎ）Ｊ　ｐｙ（Ｎ　　１　
）］／　ｈ　ｙ’　（Ｍ））＋Ｕ（Ｍ、Ｎ）　　　　　
・・・（２２−１）ｆ　ｎ（Ｍ＝Ｎ）　＝　（１／　ｑ
　）　（［Ｊ　ｎｘ（Ｍ）−Ｊ。Ｘ（Ｍ−１）］／　ｈ
　ｘ’　（Ｍ））＋（１／ｑ）　（［Ｊｎｙ（Ｎ） −Ｊ　ｎｙ（Ｎ−１）］／　ｈ　ｙ’　（Ｍ））＋Ｕ（
Ｍ、Ｎ）　　　　　・・・（２２−２）ｆ◆（Ｍ、Ｎ１
＝［１／　ｈｘ’（Ｍ）］×（［ψ（Ｍ＋１．Ｎ）−ψ
　（Ｍ、Ｎ）コ／ｈＸ（Ｍ）−［ψ（Ｍ、Ｎ）−ψ（Ｍ
−１、Ｎ）Ｍｈｘ（Ｍ−１））＋［１／ｈｙ′（Ｎ）］ ×（［ψ（Ｍ、Ｎ＋１）−ψ（Ｍ、Ｎ）コ／ｈｙ（Ｎ）
−［ψ（Ｍ、Ｎ）−ψ（Ｍ、Ｎ−１）コ／　ｈ　ｙ（Ｎ
−１））十（ｑ　／　Ｅ　）［ｒ（Ｍ、Ｎ）＋　ｐ　（
Ｍ、Ｎ）　−ｎ　（Ｍ、Ｎ）］・・・　（２２−３） 上記（ｄ）、（ｅ）、（ｆ）式を時間積分することによ
り、前記連立方程式系の解を求めるという動作解析法に
おける空間電荷の関係え□λ。、λφを安定性理論を用
いて決定することを特徴とする。

ただし、上記の説明においては、簡単のため、原関数ｆ
、、ｆｎが物理的時間τに関する微分を含まない直流定
常問題を取り上げたが、これらを含む場合については、
式（８−１）、　（８−２）において、ａｐ／ａτおよ
びａ　ｎ　／　ａτを夫々差分式％式％ の形に置き換えたものから出発すれば、本計算法の原理
に本質的な差異を生じることなく、前記連立方程式系の
解を求めることができることを指摘しておく。

本発明の第３の所定の物理現象の解析法における空間電
荷の関係λの決定法は、所定の物理現象を表わす単一の
方程式または連立方程式系ｆ（ｘ）＝Ｏをコンピュータ
により解析するために、前記方程式または方程式系を時
間微分項ｄｘ／ｄｔと空間電荷の関係λとを含む下ｖ、
（ａ４）式に置き換え、ｄｘ／ｄｔ＝λｆ（ｘ）　　　
　　　　・・・（ａ４）前記（ａ４）式を時間積分する
ことにより、前記行列方程式の解を求めることを特徴と
するる解析法において、空間電荷の関係λを安定性理論
から決定することを特徴とする。

ただし、上記において、ベクトル関数ｆは、未知量Ｘお
よびその関数の微分作用素、積分作用素、その他の演算
作用素を含んでもよいものとする。

更には、同ベクトル関数ｆが物理的時間λおよびその関
数についての微分作用素、積分作用素、その他の演算作
用素を含んでもよいものとする。

本発明の第Ｉの解析装置は、半導体デバイスモデリング
における電位と空間電荷の関係を表わすポアソン方程式
を微分方程式ｄψ／ｄｔ＝λｆφに置き換えて解析する
装置において、前記装置は、単数換えて解析する装置に
おいて、前記各論理セルに結合し調整可能な増幅利得λ
、を有するｍ個のアンプまたはこれと等価な乗算器と、
前記各論理セルに対してネットワーク状に設けられ前記
装置内の過度現象が定常状態になったことを検知する検
出手段とからなる装置において、増幅利得λ、を自動的
に決定することを特徴とする。

本発明の第２の解析装置は、半導体デバイスモデリング
における、電子、正孔の輸送方程式およびポアソン方程
式からなる連立方程式系を時間微分項と空間電荷の関係
とを含む微分方程式ｄｐ／ｄｔ＝−λｐｆ　ｐｔ　　ｄ
　ｎ　／　ｄ　ｔ　＝λｎｆｎ、　ｄψ／ｄｔ＝λｎｆ
φに置き換えて解析する装置において、前記装置は、単
数換えて解析する装置において、前記論理セルに結合し
調整可能な増幅利得λｉを有するｍ個のアンプまたはこ
れと等価な乗算器と、前記各論理セルに対してネットワ
ーク状に設けられ前記装置内の過度現象が定常状態にな
ったことを検知する検出手段とからなる装置において、
増幅利得λｉを自動的に決定することを特徴とする特本
発明の第３の解析装置は、所定の物理現象を表わす単一
の方程式または連立方程式系　ｆ（ｘ）＝Ｏの解を求め
る計算を微分方程式ｄ　ｘ　／　ｄ　ｔ　＝λｆ　（ｘ
）の積分問題に置き換えて実行するように組み立てられ
た解析装置が、時間軸上の積分プロセスを実行するアナ
ログ計算実行部と、計算途上で係数行列λの大きさを個
別に調節、適正化する手段とを有する装置において係数
行列λを自動的に決定することを特徴とする。

（作用） 本発明の方法では、半導体デバイスモデリングにおける
電位と空間電荷の関係を表わすポアソン方程式の解法、
半導体デバイスモデリングにおける電子、正孔の輸送方
程式およびポアソン方程式からなる連立方程式系の解法
などにおいて、上記各式を、時間微分項と空間電荷の関
係λとを有する微分方程式に置き換え、解く場合に、空
間電荷の関係λを自動的に設定するので、対象となる方
程式がさらに大形化しても行列解法を高速にかつ、外か
らの調整なしに実行できる。

本発明の装置は、式（ａ４）の独立変数ｔを時間とみな
し、未知量又は式（ａ４）を満足する電気回路網の所定
端子の電圧で表わされる回路網の１部分をなす。

この回路網全体の過渡現象が時間と共に進行して定常状
態が実現すれば、それがすなわち式（ａ４）の解を与え
るように装置全体が構成される。

一般に、従来法に基づくディジタル計算では。

行列方程式Ａ　ｘ　＝　ｂの解の計算において、消去法
を実行する際は、形式上ｘ　”　Ａ　−１ｂと書かれる
、所謂逆行列計算ないしは実質上これと等価な計算が含
まれ、これを実行するために多大の計算時間が消費され
る。

これに対し、本発明に用いられる方程式（ａ４）では、
逆行列の計算が含まれない。単にある現時刻の関数値ｘ
０　を用いて、行列Ａとこれとのかけ算Ａ　ｘ　’、お
よびこれから定数ベクトルｂをさし引くための引き算、
更にこの結果得られたＡｘｏ−ｂに対して、左から時間
きざみの長さδｔと係数行列λをかけたものが、回路網
の接続によりきまる時間微分量の大きさと等しくおかれ
る。この計算過程はいわば順方向の計算だけから成り立
っている。

順方向の計算が有利なことは、特に行列Ａが多数のゼロ
要素を含む場合において著るしい。この場合かけ算Ａ　
ｘ　０　の実行に際してへの非ゼロ要素のみを考慮すれ
ばよいからである。

本発明がその一部をなす。計算法が行列方程式Ａｘ＝ｂ
の真の解を与えるためには、ひとつの条件がある。それ
は、微分方程式（ａ４）の、時間軸上の積分計算が十分
長く実行された場合、未知量Ｘが定常値に到達すること
である。これを式で書けば、 Ｑｉｍ（ｄ　ｘ　／　ｄ　ｔ　）＝　Ｏｔ　−＋■ となる。この条件が満たされるものとすれば、十分な時
間が経過したあとは、微分方程式（ａ４）の左辺はゼロ
となる。従って右辺もゼロでなければならないことから
、求まったＸはＡｘ＝ｂを満足する。

上記の計算実行過程において、Ｘが定常値に到達できる
か否かをきめるのは、空間電荷の関係としての対角係数
行列λの選定の適否である。一般にλが小さ過ぎると、
時間微分量ｄｘ／ｄｔが小さすぎるため、初期試行値ｘ
＝ｘ０からの変化が殆ど起らず、従って定常解が得られ
るまでに多大の時間がかかってしまう。これに対して、
λが大きすぎると、時間微分量が大きすぎるため微分方
程式（ａ４）で表わされる方程式系（あるいは回路網）
全体の挙動が不安定になる虞がある。実際にディジタル
シミュレーションにより数値計算を行なった例をみると
、λが過大な場合、確かに未知量Ｘの変化が過大となり
、微小時間が経過した後、これを修正するめため先と逆
向きの変化がおこる。

時間の進行と共に、このことが繰り返されると、結局は
時間軸上の振動現象が起って、未知量又はいつまでも定
常値に達し得ないことが分る。未知量変化が更に大きい
と、振動の振幅が無限に大きくなるとか、または未知量
値が一方向（正または負のいずれか）に向って無限に大
きくなるなどの発散現象を生じることも知られている。

以上を要約すれば、本発明がその一部をなす計算法によ
り本発明がその一部をなす装置を適性に稼働させるため
には、対角係数行列λの大きさを適正値に決めることに
より、はどよい大きさの未知量変化を生じさせることが
重要であり、本発明はそのλの適正な決定法に関するも
のである。

更に具体的な議論をすれば、λは対角係数行列であるか
ら、ベクトルＸの成分であるｎ個の要素（Ｘｌ、　Ｘ、
、・・・ｘｎ）に対して個々にｎ個の係数（λ□、λよ
、・・・　λ。）が対応することに注意を要する。この
ため、本発明がその一部をなす計算法及び装置には、ｎ
個の未知量要素に対して個別にｎ個の係数を調節し、各
々を最適値にきめ得るという自由度が存在する。この点
において従来行列解法として知られた反復法（この方法
では所謂緩和係数が導入されるが、その個数は通常、行
列方程式全体に対して１個である）と比べて遥かに効率
よく、未知量を定常値に至らせることができる。

λと呼ばれるパラメータ関数を以下のように安定性条件
から決定する。ここでいう安定性条件とは、数値的な計
算の過程で生ずるすべての丸め誤差の累積効果が無視で
きることと定義される。（Ｇ。

Ｄ、スミス「電算機による偏微分方程式の解法」サイエ
ンス社１９７９刊のｐ６３） 先ず（ａ−４）において、変数Ｘをｘ０＋δＸで置き換
え、δＸについて線形化しδＸを含まない項を定数とし
て省略すると（ａ−５）になる。ここにδｆはｆの第一
変分である。

ｄ（δｘ）／ｄｔ＝λ（δｆ）（δｘ）−（ａ−５）（
ａ−５）を離散化すると、（ａ−６）のようになる。

δｘｔ＋ａｔ＝λｄｔ　　（δｆ）　　δＸｉ＋δｘｔ
・・・（ａ−６） このδＸｉからδＸ　ｔ＋ｃｔｔへの状態の遷移を示す
式をまとめて行列の形で書くと（ａ　−７）になりその
δＸｔの係数行列Ａを誤差伝播行列という。

δＸｔ＋ｄｔ＝Ａδｘｔ−（ａ−７） ここで関数λは必ず時間増分ｄｔとの積の形で式の中に
現れる。したがってパラメータλの決定には時間増分ｄ
ｔのきめかたが関わってくる。まず、安定性の議論に必
要な数学の定理であるＧ　ｅｒｓｃｈｇｉｒｉｎの定理
を次に示す。（Ｇ、Ｄ、スミス「電算機による偏微分方
程式の解法」サイエンス社１９７９刊のｐ６８） 定理（Ｇｅｒｓｃｈｇｉｒｉｎ　ｃｉｒｃｌｅ　ｔｈｅ
ｏｒｅ＋ｍ）　　：ｎ次元行列Ａ　＝　（ａ工、）の任
意の固有値μは、Ａ＝Ｇ＋Ｆ、　Ｆ＝（４１ｊ）、　Ｇ
＝ｄｉａｇ（ａ１□、　ａ、□。

・・・ｙ　ａ　ｎｎ）とすると、次のような閉円板の集
合Ｇ。

に含まれる。

すなわち μＣｕ　Ｇ。

ただし　Ｇ４＝（ｚは複素数；ＩＺ−ａｉｉｌ≦ｒｔ）
−（ａ−８）ｒ工＝Σ１ｆｉＮ この定理は複素平面上での行列の固有値の分布を規定す
るものである。これを（ａ−７）に適用する。安定性の
条件、つまり固有値の絶対値が１以下であるという条件
は、（ａ−８）で示される円板群が、複素平面上の単位
円内に存在すれば満足される。この条件を満たすように
λを決定する式を陽な形で書下だし、本発明がその一部
を構成する回路網の装置を組み立てるに際しては、この
λに対応する増幅利得を持つようにする。λを陽な形で
書下せない場合には、数値的又は記号的にλを決定する
ようにしても良い。

実際に本発明の装置を組立てるに際しては、後述のよう
にｎ個の論理セルの各々に付随したアンプを設け、その
各々が増幅利得λｉをもつよう構成する。各アンプは個
々の実際問題に十分対応可能な程度の十分大きな増幅利
得をもつように設計し、利得調整用のつまみまたはコン
トローラを具えるものとする。

（実施例） 本発明の第１実施例に係る半導体デバイスの動作解析法
における空間電荷の関係の決定法を説明する。

デバイスモデリングにおいて、例えばｐｎダイオードの
ｐ側をゼロ電位、ｎ側を正電位■に保つ場合、ｐｎ接合
の界面近傍に空乏層が形成される。

この際、電位■がダイオードの耐圧以下であれば、ダイ
オードに流れる電流はほとんど無視できる程度に小さい
。このような場合、基本方程式（１）〜０のうちポアソ
ン方程式■だけを解けばデバイス全体の電位φの分布が
正確に求まることが知られている（例えば、Ｋｕｒａｔ
ａ、　”Ｎｕｍｅｒｉｃａｌ　Ａｎａｌｙｓｉｓｆｏｒ
　Ｓｅｍ１ｃｏｎｄｕｃｔｏｒ　Ｄｅｖｉｃｅｓ、　Ｈ
ｅａｔｈ、　１９８２゜第８章、又は、倉口、「バイポ
ーラトランジスタの動作理論」、近代科学社、昭和５５
年３．２節および５．１節参照）。この際正孔密度およ
び電子密度Ｐ。

ｎは各々次式で与えられる。

ｐ＝ｎ４ｅ−θφ　　　　　　　　　・・・（２４−１
）。＝ｎ□ｅ＄（＋−Ｖ）　　　　　　　　・・・（２
４−２）これらを弐〇に代入すれば、２次元空間の場合
法の方程式が得られる。

Ｃａ”ψ／ａｘ２）＋（ａ２ψ／ａｙ”）＝−（ｑ／ε
）　［ｒ（ｘ＋　ｙ）　＋　ｎｔｅ−“−ｎｉ　ｃｅ（
＄−Ｖ）コ　　　　　　　　　　　　　　　　　　　　
・・・　（２５）この方程式を離数系に変換することは
、先と同じ要領で行われる。その具体的な形は、式■の
右辺に（２４−１）、　（２４−２）を代入したものに
他ならない。すなわち、 ［１／　ｈ　ｘ’　（Ｍ）］ ×（［ψ（Ｎ＋１．Ｎ）−ψ（Ｍ、Ｎ）］／ｈｘ（Ｍ）
−［ψ（Ｍ、Ｎ）−ψ（Ｍ−１，Ｎ）コ／　ｈ　ｘ（Ｍ
　　　１　））＋　［１／　ｈ　ｙ’　（Ｍ）］ ×（［ψ（Ｍ、Ｎ＋１）−ψ（Ｍ、Ｎ）コ／ｈｙ（Ｎ）
−［ψ（Ｍ、Ｎ）−ψ（Ｍ、Ｎ−１）コ／ｈｙ（Ｎ−１
））＝−ρ（Ｍ、Ｎ）／ε　　　　　　　　　　　　・
・・（２６）ただし ｐ　（Ｍ＋Ｎ）＝　ｑ　［ｒ（Ｍ、Ｎ）＋　ｎ　Ｌ　ｅ
　−”””’−□、　ｅθ（φ（Ｍ、Ｎ）−Ｖ）コ （２７） である。従来法によりこの方程式の解を求めることは、
前に詳述したものと同じ要領に従い、行列ベクトル方程
式の解の計算と反復計算の実行に帰着する。

これに対し本発明をその一部として含む方法では、時間
微分項と空間電荷の関係λを含むつぎの形の方程式から
出発する。

ａ　ψ（Ｍ、Ｎ）／ａ　ｔ＝ｔ（Ｍ、Ｎ）ｘ　（［１／
ｈｘ’（Ｍ）］ Ｘ（［ψ（Ｎ＋１．Ｎ）−ψ（Ｍ、Ｎ）］／ｈｘ（Ｍ’
）−［ψ（Ｍ、Ｎ）−ψ（Ｍ−１，Ｎ）／ｈＸ（Ｍ’−
１））＋　［１／　ｈ　ｙ’　（Ｎ）コ ×（［ψ（Ｍ、Ｎ＋１）−ψ（Ｍ、Ｎ）］／ｈｙ（Ｎ’
）−［ψ（Ｍ、Ｎ）−ψ（Ｍ、Ｎ−１）］／ｈｙ（Ｍ’
−１））＋ρ（Ｍ、Ｎ）／Ｆ）　　　　　　　　　　　
　　・・・（２８）この方程式を時間軸上に沿って積分
する。

そのためのフローチャートを第３図に示す。即ち、本方
法をその１部として含む方法においては、まず計算をス
タートするにあたって、ステップａ２において、半導体
デバイスの構造、不純物データをコンピュータのメモリ
に入力する。次にステップｂ２において、半導体デバイ
スの分割点配置を決定すると共に、ステップｃ２でバイ
アス電圧■を入力する。次に電位の初期試行値ψ０をス
テップｄ２で与え、ステップｅ２において時間積分開始
のための初期時間ｔ１を定める。

次にステップｆ２において、ψ１を計算する。

ψ１は次式で与えられる。

ψ１（Ｍ、Ｎ）＝ψ’（Ｍ、Ｎ） ＋δｔ・λ（Ｍ、Ｎ） Ｘ　（［１／ｈＸ’（Ｍ）］ ×（［ψ’（Ｎ＋１．Ｎ）−ψ’（Ｍ、Ｎ）］／ｈＸ（
Ｍ’）−［ψ’（Ｍ、Ｎ）−ψ’（Ｍ−１，Ｎ）／ｈＸ
（Ｍ’−１））＋［ｌ／ｂｙ’（Ｎ）コ Ｘ（［ψ’（Ｍ、Ｎ＋１）−φ’　（Ｍ、Ｎ）］／　ｈ
　ｙＯｔｙ’）−［ψ’（Ｍ、Ｎ）−ψ０（Ｍ、Ｎ−１
）］／ｈｙ（Ｍ′−１））＋ρ（Ｍ、Ｎ）／ε）　　　
　　　　　・・・（２９）式（２９）に従って右辺の値
をすべての点１≦Ｍ≦Ｋ、１≦Ｎ≦Ｌについて求め、こ
れらをψ１（Ｍ、Ｎ）；ｌ≦Ｍ≦に、１≦Ｎ≦Ｌに等し
いとおく。これによって時間軸上の積分が１ステップ進
行したことになる。

次にステップｇ２において、変量の変化分の絶対値１ψ
１（Ｍ、Ｎ）−ψ’（Ｍ、Ｎ）ｌを全点について求め、
これらの全てが所定の誤差限界より小さいか否かをチエ
ツクする。

もしデバイス空間の中に１点でも誤差限界より大きな変
化を生じた点があれば、ステップｈ２において、先に求
めた修正値ψ’（Ｍ、Ｎ）を新な試行値とみなし、これ
らをψ′５（Ｍ、Ｎ）に等しいとおき、更に時刻ｔ１を
ｔｏとして、前述と同じ要領で更に時間軸上の積分計算
を行う。

もしすべての点で変量変化の絶対値が誤差限界を下まわ
ったならば、定常解が得られたものとみなして計算を終
了する。その結果得られた解は、（２１−３）において
、式の左辺の時間微分項がゼロになるようなものである
から、ｆ◆（Ｍ、Ｎ）＝Ｏ；１≦Ｍ≦に、１≦Ｎ≦Ｌが
威り立つ。このことは即ち、方程式■あるいは方程式（
２６）の解が求まったことに他ならない。

なお、上記において時間軸上の積分を実行する際に、未
知変量ψ（Ｍ、Ｎ）の変化の絶対値１ψ１（Ｍ、Ｎ）−
ψ’（Ｍ、Ｎ）ｌが小さすぎると、計算過程の安定性は
よいものの、時間ｔが経過してもψ（Ｍ、Ｎ）はなかな
か真の解に到達できない。この反面、もし上記の変化の
絶対値が大きすぎると。

計算過程が不安定となり、時間ｔの経過に伴ってψ（Ｍ
、Ｎ）の振動または発散を生じるから、真の解が求まら
ない。

本発明で決定法を示す空間電荷の関係λ（Ｍ、Ｎ）は。

このような変量変化を適度に生じさせるための調整因子
である。特に１本例のように半導体デバイスのポアソン
方程式で、かつ式（２９）に示した２次元空間モデルの
場合、上記の条件を満たすようなλ（Ｍ、Ｎ）を次のよ
うに定める。先ず（２８）において、変数ψをψ０＋δ
ψ０で置き換えたのちδψについて線形化する。そして
、線形化した式について安定性を考える。ここでδψを
含まない項を定数として省略する。すると、 δψ１（Ｍ、Ｎ）＝δｔ・λ（Ｍ、Ｎ）（Ａδψ’（Ｍ
、Ｎ−１）十Ｂδψ’（Ｍ−１、Ｎ）＋Ｃδψ’（Ｍ、
　Ｎ＋１）＋Ｄδψ’（Ｍ＋　１　、Ｎ）−Ｈδψ”（
Ｍ、Ｎ）−αδψ’（Ｍ、Ｎ））＋δψ’（Ｍ、Ｎ）（
３０−１） Ｈ＝１／［ｈｘ（Ｍ’　　１）ｈｘ’（Ｍ）］］＋／［
ｈｘ（Ｍ’）ｈｘ’（Ｍ）］ ＋　１／［ｈｙ（Ｎ’　　　１）ｈｙ’（Ｎ）コ＋　１
／［ｈｙ（Ｎ’）ｈｙ’（Ｎ）コα＝（ｑθｎｌ／ε）
× ［ｅｘｐ（−θψ’（Ｍ、Ｎ）） ＋ｅｘｐ（θ（ψ０（Ｍ、　　Ｎ）−Ｖ））コ・・・　
（３０−６） （３０−７） 以上をまとめてベクトルの形で書くと ψ１＝Ａψ　　　　　　　　　　・・・（３０−８）に
なる。

ここでＧ　ｅｒｓｃｈｇｏｒｉｎの定理を利用すると、
安定性の条件、つまり誤差伝播行列への固有値の絶対値
が１以下であるという条件は、Ｇ　ｅｒｓｃｈｇｏｒｉ
ｎの定理で示される円板群の中心座標は実数であるから
、円板群の中心は実軸上にあるため、実軸でのみ考えて
、次の２つの式が満たされれば良い。

１−λ（Ｍ、Ｎ）・δｔ（Ｈ＋α） −λ（Ｍ、Ｎ）・δｔ−Ｈ＞−１・・・（３０−９）↓
−λ（Ｍ、Ｎ）・δｔ（Ｈ十α） ＋λ（Ｍ、Ｎ）・δｔ−Ｈ＜１　　　・・・（３０−１
０）上記のうち（３０−１，０）は−λ（Ｍ、Ｎ）・δ
ｔ・αくＯとなり自明であるから、（３０−９）の不等
式が安定性条件を与える。すなわち（２９）が安定であ
るためには、時間増分Δｔと関数λを、空間刻みＨと比
線形項の大きさαの関数として次式によりきめればよい
。

ここで、空間刻みＨは、場所のみの関数で時間には依存
しないのに対し、非線形項αはψを含むから場所・時間
の関数である。従ってλは結果的に場所と時間の関数に
なる。

以上の議論を１次元、２次元に適用した場合、変化する
のはＨの表示式のみである。また、差分法以外の離散化
方法（有限要素法、有限体積法など）でも、以上のλの
決定法の議論はそのまま適用できる。

実際に利用する場合には、式（３０−１１）の右辺の値
より少し小さめにλ・δｔをきめるので、そのためのパ
ラメータω（ωは１以下の１に近い正数）を利用し以下
のように定義する。

以下、式（３０−１２）の意味を定性的に説明する。

一般に半導体デバイスの内部で不純物濃度が大きい領域
では電子、正孔密度ｎｅＰのいずれか一方が不純物濃度
とほぼ相等しく、空間電荷中性条件が保たれる。この場
合αは、Ｈおよび分子の２に比べて非常に大きいので、
λ（Ｍ、Ｎ）・δｔは非常に小さい。の結果、変量変化
１ψ１（Ｍ、Ｎ）−ψ’（Ｍ、Ｎ）ｌ　　もまた非常に
小さい。これに対して、不純物濃度が比較的小さく、空
乏層が形成されやすい領域では、αはＨと同程度か、ま
たはＨに比べて無視できる程度に小さい。その結果λ（
Ｍ、Ｎ）・δｔは比較的大きく、従って、変量変化も比
較的大きい。

以上の理由により、計算開始のための初期条件として、
変量ψ（Ｍ、Ｎ）に対し、デバイス中の至る所で空間電
荷中性条件を満足するような値を与えれば、時間ｔの経
過と共に低濃度領域では大きな変量変化が生じて、空乏
層の形成が促進され、変量ψ（Ｍ、Ｎ）は真の解に向か
って収束することになる。

この計算を通常の如く大形計算機を用いたディジタル計
算により実行することは勿論可能であるが、他に、電子
回路を応用することにより、その回路の時間応答（過渡
現象）が定常状態に達したならば、その状態がすなわち
所定の方程式の解に対応するような装置を構築すること
が可能である。

この点については以下に本発明の装置として説明する。

第１の実施例の方法を実行する装置について説明する。

このような装置を組み立てるには、Ｎ丁個の論理セルの
各々に付随したアンプが増幅利得λｉを有するように構
成する。ここでＮ丁は、連続系に属する原方程式をコン
ピュータを用いて解くために、差分法、有限要素法など
の近似解法に基づき、原方程式を離散系に置き換えた場
合のメツシュの分割点の総数である。各アンプは個々の
実際問題に十分対応可能な程度の大きな増幅利得を有す
るように設計する。このとき、利得は（３〇−１２）に
示されるように、入力量から加算器／乗算器／除算器を
用いて決定されるものとする。全回路網で１つ決定され
ている（３０−１２）のパラメータωを調整するつまみ
又はコントローラを備えるとする。

ここで利得を計算する式（３０−１２）は、α又はＨの
値によって区分的に近似式を用いることも含める。また
、利得は入力量から回路をもって決定されるものを含む
。

ただし、問題により増幅利得は１より小さいこともあり
得るので、この際には、アンプは実際上減衰器となるが
、このような場合も包括されるものとする。

半導体デバイスモデリングにおける、電位と空間電荷の
関係を表すポアソン方程式をコンピュータにより解析す
る方法について具体的に実施した数値計算の結果を示す
。計算対象は第４図に示すような不純物分布・寸法をも
つｐｎ接合を持つデバイスである。これを３次元問題と
して解析する場合は２０　ｘ　２６　ｘ　４２メツシユ
（２１８４０自由度）、２次元解析の場合そのｘ−ｚ平
面による切断面を２０ｘメツシユ（８４０自由度）を採
用する。ゼロバイアスでの計算の試行値は、次のように
与えている。

ここで本発明をその一部として含む方法において、本発
明が述べたような方法で空間電荷の関係λを与えた時の
陽的反復の回数と収束の状態（３次元）を表１と第５図
（ａ）〜（ｅ）に示す。ここで誤差Ｗは、 Ｗ＝關δψ／（ψ０＋δψ）Ｉｌ− とし、これが１．Ｅ−８になるまで反復を繰り返すとす
る。結局、 １）ωが小さすぎると収束に要する回数が多くなる。ω
が大きくなるにつれ収束は早くなるが。

太きすぎる（ω＞０．９　”）と急激に逆に収束しにく
くなる。

２）ω＜１．０の場合には反復の初期にはＷに大きい変
動が見られるものの、比較的早期にこのパターンから脱
出して収束はほぼ単調になる。

３）ω＝１．０の場合にはＷが小さくなる方向であるが
、収束にいたらない。

４）ω〉１．０では収束に向かわない。またＷの値に不
規則な振動がある。

が分かる。この場合、収束が最も速いのはω：０．９の
ときである。

表工　収束に必要な反復回数（３次元、ゼロバイアス）
計算対象が２次元、ゼロバイアスの場合の結果は表２．
第５図（ｇ）のようになる。収束の特性等については３
次元、ゼロバイアスの場合と同様である。この場合、最
も収束の早いωの値は、ω＝０．８８であり３次元の場
合と若干ずれている。この表から、腸内反復では、ωを
どの様に選ぶかが計算時間に効くが、はぼω＝０．９程
度にすればωがこれにより小さめにずれた場合でも収束
回数の変化は小さいから、事実上問題ない。

表２　収束に必要な反復回数（２次元、ゼロバイアス）
次に、従来の陰解法を利用した方法と、本発明をその一
部として含む方法において、本発明が述べたような方法
で空間電荷の関係を与えた時との計算時間の比較につい
て述べる。比較は２次元、３次元のゼロバイアスの時の
計算を使用した。従来の陰解法ではニュートン法により
非線形方程式を解いており、その際大型線形連立方程式
解法としてはガウスの消去法／ＢＣＧ法ｌＣＲ法を利用
している。表３に結果を示す。これより２次元、３次元
の計算に於いて、本発明をその一部として含む方法は従
来法に比べて２〜４倍高速であることがわかる。この比
較を計算精度的に見ると、ニュートン法による方法は初
め誤差が大きくても収束し始めると急速に収束する。そ
れに対して腸内反復では、比較的緩慢に解に近付いてい
く。従って収束判定の許容誤差Ｗを大きくする（例えば
１．Ｅ−４）と更に腸内反復が有利になる。

表３　計算時間の比較 次に本発明をその一部として含む方法による逆バイアス
計算をした結果について述べる。図でｖＲを１ｖから２
００Ｖまで順次上げた。電圧をあげた時の計算の試行値
は、その−段階前のバイアス電圧ｖ１での解とかける電
圧の差Ｖ２−Ｖｌとをもとにして次ように決めている。

第６図に逆バイアス計算での残差と腸内反復回数の関係
を示した。電圧が上がると収束に要する反復回数が増加
するものの、いずれの場合も収束し解が得られた。これ
らの計算結果から逆バイアス計算にも本発明をその一部
として含む方法が利用できることが示された。

本発明の第２実施例に係る半導体デバイスの動作解析法
における調整因子の決定法について述べる。

本発明では、その１部として含む方法では方程式（８−
１）、　（８−２）、■）の解を求めるにあたり、これ
らを時間微分項を含む、つぎの形におきかえる。

ｄ　ｐ　（Ｍ、Ｎ）／　ｄ　ｔ　＝−λｐ（Ｍ＝Ｎ）ｆ
ｐ（Ｍ、Ｎ）・・・　（２１−１） ｄ　ｎ　（Ｍ、Ｎ）／ｄ　ｔ　＝λｎ（Ｍ、Ｎ）ｆｎ（
Ｍ、Ｎ）・・・（２１−２） ｄψ（Ｍ、Ｎ）／ｄ　ｔ　＝−λφ（Ｍ、Ｎ）ｆ◆（Ｍ
、Ｎ）・・・　（２１−３） ただし各式の右辺のｆＰ、ｆｎ、ｆφは次式のとおり定
義するものとする。

ｆｐ（Ｍ、　Ｎ）＝（１／ｑ）　（［ＪＰＸ（Ｍ）Ｊ　
ｐｘ（Ｍ　　１　）］／　ｈ　ｘ’　（Ｍ））＋（１／
ｑ）　（［Ｊｐｙ（Ｎ） Ｊ　ｐｙ（Ｎ−１）］／　ｈ　ｙ’　（Ｍ））＋Ｕ（Ｍ
、　Ｎ）　　　　・・・（２２−１）ｆｎ（Ｍ、　Ｎ）
＝（１／（Ｉ）　（［Ｊｎｘ（Ｍ）−Ｊ　ｎｘ（Ｍ−１
）］／　ｈ　ｘ’　（Ｍ）十〇　／　ｑ　）　（［Ｊ　
ｎｙ（Ｎ）Ｊ　ｎｙ（Ｎ−１）］／　ｂ、　ｙ’　（Ｍ
））＋Ｕ（Ｍ、　Ｎ）　　　　・・・（２２−２）ｆφ
（Ｍ、Ｎ）　　＝（［１／ｈＸ’（Ｍ）］×（［ψ（Ｍ
、Ｎ）−ψ（Ｍ、Ｎ）コ／ｈＸ（Ｍ）−［ψ（Ｍ、Ｎ）
−ψ（Ｍ−１，Ｎ）コ／　ｈ　ｘ（Ｍ−１））＋［１／
ｈｙ’（Ｎ）］ ×（［ψ（Ｍ、Ｎ＋１）−ψ（Ｍ、Ｎ）］／ｈ　ｙ（Ｎ
）−［ψ（Ｍ、Ｎ）−ψ（Ｍ、Ｎ−１）］／ｈｙ（Ｎ−
１））＋（ｑ／ε）［ｒ（Ｍ、Ｎ）＋ｐ（Ｍ、Ｎ）−ｎ
　（Ｍ、　Ｎ）］　　・・・（２２−３）再び（２１−
１）〜（２１−３）にもどり、各式の左辺の時間微分項
を差分近似する。この際、時間軸上に有限個の分割点を
設け、旧時刻の値は上つきの数字Ｏ１１部刻の値は上つ
きの数字工で表わすことにする。その結果つぎの各式が
得られる。

ｐ’（Ｍｅ　Ｎ）＝ｐ’（Ｍｅ　Ｎ） −δｔλＰＣＭ、　Ｎ）ｆＰ。（Ｍ、　Ｎ）（２３−１
） ｎｌ（Ｍ、Ｎ）＝ｎ’（Ｍ、Ｎ） 一δｔλｎ（Ｍ、Ｎ）ｆ、’（Ｍ、Ｎ）（２３−２） ψ１（Ｍ、Ｎ）＝ψ’（Ｍ、Ｎ） −８ｔ　λ、（Ｍ、Ｎ）ｆ、ａ（Ｍ、Ｎ）・・・（２３
−３） ただし、１部Ｍ≦Ｋ、工≦Ｎ≦Ｌとする。

以下、第７図のフローチャートを参照して説明する。本
方法をその１部として含む方法においては、まず計算を
スタートするにあたって、ステップａｌにおいて、半導
体デバイスの構造、不純物データをコンピュータのメモ
リに入力する。次にステップｂ１において、半導体デバ
イスの分割点位置を決定すると共に、ステップＣ１でバ
イアス電圧Ｖを入力する。次に基本変量ｐ（Ｍ、Ｎ）、
ｎ（Ｍ、Ｎ）、ψ（Ｍ、Ｎ）の初期試行値ｐ’（Ｍ、Ｎ
）。

ｎ　’　（Ｍ　、　Ｎ　）　、ψ’（Ｍ、Ｎ）をステッ
プｄ　ｉで与え、ステップｅ１において時間積分開始の
ための初期時間ｔ１を定める。

次にステップｆ１において、式（２３−１）〜（２３−
３）に従って各式の右辺の値をすべての点１≦Ｍ≦Ｋ、
１≦Ｎ≦Ｌについて求め、これらをｐ’（ＭｙＮ）−ｎ
ｌ（ＭｔＮＬ　ψ１（Ｍ、Ｎ）；　１部Ｍ≦Ｋ、１≦Ｎ
≦Ｌに等しいとおく。これによって時間軸上の積分が１
ステップ進行したことになる。

次にステップｇ１において、各変量の変化分の絶対値ｌ
　ｐ”（ＭｔＮ）−ｐ’（ＭｔＮ）　ｌ　ｔ　Ｉ　ｎｌ
（Ｍ−Ｎ）−ｎ’（Ｍ、Ｎ）　ｌ　、　ｌψ”（Ｍ、Ｎ
）−ψ’（Ｍ、Ｎ）ｌを全点について求め、これらの全
てが所定の誤差限界により小さいか否かをチエツクする
。

もしデバイス空間の中に１点でも誤差限界より大きな変
化を生じた点があれば、ステップｈ１において、先に求
めた修正値ｐ　’　（Ｍ　ｒ　Ｎ　）　＋　ｎ　１（Ｍ
　？Ｎ）、ψ１（Ｍ、Ｎ）を新な試行値とみなし、これ
らをｐ’（Ｍ、Ｎ）、ｎ’（Ｍ、Ｎ）、ψ’（Ｍ、Ｎ）
に等しいとおき、更に時刻ｔ１をｔｏとして、前述と同
じ要領で更に時間軸上の積分計算を行う。

もしすべての点で各変量変化の絶対値が誤差限界を下ま
わったならば、定常解が得られたものとみなして計算を
終了する。その結果得られた解は、（２１−１）〜（２
１−３）において、各式の左辺の時間微分項がゼロにな
るようなものであるから、ｆ　ｐ（Ｍ　、　Ｎ）”　Ｏ
、ｆ　ｎ（Ｍ、　Ｎ）＝　０　、　ｆ　＋ＣＭ、　Ｎ）
＝Ｏ；１≦Ｍ≦に、１≦Ｎ≦Ｌが成り立つ。このことは
即ち、方程式（８−ＩＬ　（ｓ−２）および■）の解が
求まったことに他ならない。

この方法における空間電荷の関係λＰ（Ｍ、Ｎ）、λｎ
（Ｍ。

Ｎ）、λφ（Ｍ、Ｎ）の決定方法は、すでに第１の実施
例で述べたことと、同じように出来る。つまり線形化し
、下記（２３−４）になり、その誤差伝播行列への固有
値が複素平面上の単位円の中にすべて含まれるという条
件からλ□λ。、λφを決定する。

〔発明の効果〕

以上に示した空間電荷の関係または増幅利得を決定する
方法を利用することで、この空間電荷の関係または増幅
利得をその一部として含む解析方法および解析装置が効
率良く実現できる。そして、この空間電荷の関係または
増幅利得その一部として含む解析方法、および解析装置
は、従来の大型の行列方程式を解いて行く方法に比べて
非常に効率が良い。

【図面の簡単な説明】

第１図は半導体デバイス空間のメツシュ分割を示す図、
第２図は電子・正孔の輸送方程式およびポアソン方程式
からなる連立方程式系を離散化したものを行列・ベクト
ルで示した図、第３図は本発明の第一実施例に関わる半
導体デバイスの動作解析法を示すフローチャートを示す
図、第４図は、本発明の第一実施例に関わる半導体デバ
イスの１つ目の計算実例の解析対象を示す図、第５図は
本発明の第一実施例に関わる半導体デバイスの１つ目の
計算実例の収束挙動を示した図、第６図は本発明の第一
実施例に関わる半導体デバイスの計算実例２つ目の計算
実例である逆バイアス計算の収束挙動を示す図、第７図
は第二実施例に関わる半導体デバイスの動作解析法を示
すフローチャートを示した図である。

Claims (6)

【特許請求の範囲】
1. （１）半導体デバイスモデリングにおける、電位と空間
   電荷の関係を表わすポアソン方程式をコンピュータによ
   り解析する方法において、前記ポアソン方程式を時間微
   分項∂ψ／∂ｔと感度係数λを含む下記（ａ）式に置き
   換え、 ∂ψ／∂ｔ＝λ＿ψｆ＿ψ・・・（ａ） 半導体デバイスの分割点配置（Ｍ、Ｎ）を決定すると共
   に、感度係数λにも前記半導体デバイスの物理的特性に
   合致した適正な空間位置依存性を与えて、上記（ａ）式
   を下記の（ｂ）式に変換し、∂ψ（Ｍ、Ｎ）／∂ｔ＝λ
   （Ｍ、Ｎ）×ｆ＿ψ（Ｍ、Ｎ）・・・（ｂ）上記（ｂ）
   式を線形化し、時間について離散化し、且つ空間位置依
   存性をまとめてベクトル表示した下記式（ｃ）式∂ψ＿
   ｔ＿＋＿ｄ＿ｔ＝Ａ∂ψ＿ｔ・・・（ｃ）における誤差
   伝播行列への固有値が１より小さくなるように感度係数
   を決めて時間積分することにより、前記ポアソン方程式
   の解を求めることを特徴とする半導体デバイスの動作解
   析法。
2. （２）半導体デバイスモデリングにおける、電子、正孔
   の輸送方程式およびポアソン方程式からなる連立方程式
   系をコンピュータにより解析する方法において、前記連
   立方程式系を、時間微分項ｄｐ／ｄｔ、ｄｎ／ｄｔ、ｄ
   ψ／ｄｔと感度係数λ＿ｐ、λ＿ｎ、λ＿ψとを含む下
   記（ａ）、（ｂ）、（ｃ）式に置き換え、 ｄｐ／ｄｔ＝λ＿ｐｆ＿ｐ・・・（ａ） ｄｎ／ｄｔ＝λ＿ｎｆ＿ｎ・・・（ｂ） ｄψ／ｄｔ＝λ＿ψｆ＿ψ・・・（ｃ） 半導体デバイスの分割点配置（Ｍ、Ｎ）を決定すると共
   に、感度係数λ＿ｐ、λ＿ｎ、λ＿ψにも前記半導体デ
   バイスの物理的特性に合致した適正な空間位置依存性を
   与えて、上記（ａ）、（ｂ）、（ｃ）式を下記の（ｄ）
   、（ｅ）、（ｆ）式に変換し、 ｄｐ（Ｍ、Ｎ）／ｄｔ＝−λ＿ｐ（Ｍ、Ｎ）×ｆ＿ｐ（
   Ｍ、Ｎ）・・・（ｂ）ｄｎ（Ｍ、Ｎ）／ｄｔ＝λ＿ｎ（
   Ｍ、Ｎ）×ｆ＿ｎ（Ｍ、Ｎ）・・・（ｅ）ｄψ（Ｍ、Ｎ
   ）／ｄｔ＝λ＿ψ（Ｍ、Ｎ）×ｆ＿ψ（Ｍ、Ｎ）・・・
   （ｆ）上記（ｄ）、（ｅ）、（ｆ）式を線形化し、時間
   について離散化し、且つ空間位置依存性をまとめてベク
   トル表示した下記（ｇ）式 ▲数式、化学式、表等があります▼…（ｇ） における誤差伝播行列Ａの固有値が１より小さくなるよ
   うに感度係数を決めて時間積分することにより、前記連
   立方程式系の解を求めることを特徴とする半導体デバイ
   スの動作解析法。
3. （３）所定の物理現象を表わす単一の方程式または連立
   方程式ｆ（ｘ）＝０をコンピュータにより解析するに際
   し、前記方程式を時間微分項 ｄｘ／ｄｔと感度係数λとを含む下記（ａ）式に置き換
   え、 ｄｘ／ｄｔ＝λｆ（ｘ）…（ａ） 感度係数λに前記物理現象に適した値を与えて、前記（
   ａ）式を線形化し、時間及び物理量について離散化しベ
   クトル表示した下記（ｂ）式 ｄｘ＿ｔ＿＋＿ｄ＿ｔ＝Ａｄｘ＿ｔ…（ｂ）における誤
   差伝播行列Ａの固有値が１より小さくなるように感度係
   数を決めて時間積分することにより前記方程式の解を求
   めることを特徴とする所定の物理現象の解析法。
4. （４）半導体デバイスモデリングにおける、電位と空間
   電荷の関係を表わすポアソン方程式を時間微分項と感度
   係数を含む微分方程式δψ／δｔ＝λ＿ψｆ＿ψに置き
   換えて解析する装置において、前記装置は、単数または
   複数のｍ個の論理セルと、前記各論理セルに結合し調整
   可能な増幅利得λ＿ｉを有するｍ個のアンプまたはこれ
   と等価な乗算器と、前記各論理セルに対してネットワー
   ク状に設けられ前記装置内の過度現象が定常状態になっ
   たことを検知する検出手段とからなり、且つ過度現象に
   おける特定時点での誤差の誤差伝播行列の固有値が１よ
   り小さくなるように増幅利得λ＿ｉを決めることを特徴
   とする半導体モデリングに使用する装置。
5. （５）半導体デバイスモデリングにおける、電子、正孔
   の輸送方程式およびポアソン方程式からなる連立方程式
   系を時間微分項と感度係数とを含む下記微分方程式ｄｐ
   ／ｄｔ＝−λ＿ｐｆ＿ｐ、ｄｎ／ｄｔ＝λ＿ｎｆ＿ｎ、
   ｄψ／ｄｔ＝λψｆψに置き換えて解析する装置におい
   て、前記装置は、単数または複数のｍ個の論理セルと、
   前記各論理セルに結合し調整可能な増幅利得λ＿ｉを有
   するｍ個のアンプまたはこれと等価な乗算器と、前記各
   論理セルに対してネットワーク状に設けられ前記装置内
   の過度現象が定常状態になったことを検知する検出手段
   とからなりかつ、過度現象における特定時点での誤差の
   誤差伝播行列の固有値が１より小さくなるように増幅利
   得λ＿ｉをきめることを特徴とする半導体モデリングに
   使用する装置。
6. （６）所定の物理現象を表わす単一の方程式または連立
   方程式系ｆ（ｘ）＝０の解を求める計算を微分方程式ｄ
   ｘ／ｄｔ＝λｆ（ｘ）の積分問題に置き換えて実行する
   ように組み立てられた解析装置が、時間軸上の積分プロ
   セスを実行するアナログ計算実行部と、計算途上で係数
   行列λの大きさを個別に調節、適正化する手段とを有し
   かつ、過度現象における特定時点での誤差の誤差伝播行
   列の固有値が１より小さくなるように感度係数をきめる
   ことを特徴とする所定の物理現象のモデリングに使用す
   る装置。