JPH03179761A - 半導体デバイスの動作解析法およびそのために使用する装置 - Google Patents

Abstract

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

Description

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

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

ａｐ／ａ　ｒ−（−１／ｑ）（ｄｌｖ　Ｊ、）＋Ｇ−、
Ｕ　　　　　　　　　　・・・　（１）９　ｎ／ａ　ｒ
−（１／ｑ）　　（ｄｌｖ　Ｊ　、　）　＋Ｇ−Ｕ・・
・　（２） Ｊｊ　−−ＱＤｐ　　（ｇｒａｄｐ）　−Ｑμｓ　ｐ　
（ｇｒａｄψ）・・・　（３） Ｊａ　　＝ｑＤ、　　（ｇｒａｄｎ）　　−４μｍ　　
ｎ　　（ｇｒａｄψ）・・・　（４） ｄｉｖ（ｇｒａｄ　ψ　）　　−（−ｑ／　　ε　）　
　　ＣＮａ　　　　　Ｎ、　　＋　ｐ−ｎ）　　　　　
　　　　・・・　（５）これらの方程式は、従来法にお
いても、本発明の方法においても共通に用いられるもの
である。

式（１）および（２）は各々正孔、電子の連続方程式を
表わす。Ｇ、Ｕは各々キャリアの発生、消滅を表わす。

式（３）、（４）は各々正孔電流密度、電子電流密度を
具体的に記述したものであり、正孔電流密度Ｊ、および
電子電流密度Ｊ７が、それぞれキャリア密度の傾斜ｇｒ
ａｄ　ｐ　、　ｇｒａｄ　ｎに比例する拡散電流成分（
第１項）と、電位傾斜ｇｒａｄψと各キャリア密度に比
例するドリフト成分（第２項）の代数和からなることを
表わす。式（５）は、電位ψと空間電荷密度の関係を表
わすポアソン方程式である。

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

方程式（１）〜（５）の数値解を求めるための第１のス
テップは、これらの式を連続形（微分方程式）から離散
系（差分方程式）に書き換えることである。この際、解
析の対象となるデバイス空間を長方形のメツシュに分割
する必要がある。二次元空間モデルについてのメツシュ
分割を概念的に第９図に示す。縦、横の分割点間隔り、
、ｈ、を第９図のように定義し、かつ、各々の中間点の
距離り、、ｈ、　　を次式で定義する。

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

ｄｌｖ　Ｊ、　−ａａ　Ｊａｍ／ａ　ｘ＋Ｅ）　Ｊｐｙ
／ａ　ｙ＝［Ｊ　　□（Ｍ’）−Ｊ　　、、（Ｍ’−１
）コ／ｈ、　　　　（Ｍ） ＋　　［Ｊ　、、（Ｎ’）−Ｊ　、、（ｆ−ｔ）］／ｈ
、’（Ｎ） ・・・　（７） 電子電流についても同様の書き替えが可能であるが、自
明であるので省略する。

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

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

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

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

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

ここで、以下の展開のために、キャリア再結合項Ｕ　（
Ｍ、Ｎ）を具体的に５ｈｏｃｋ　ｌ　ｅｙ−Ｒｅａｄ−
Ｈａ　ｌ　ｌの形で与える。即ち、 Ｕ　（Ｍ、Ｎ）＝　［ｐ　（Ｍ、Ｎ）　　・ｎ　（Ｍ、
Ｎ）ｎ　　＋　　２　　（Ｍ、　　　Ｎ）　　コ　／　
　（ｒ、　　　［ｐ　　（Ｍ、　　　Ｎ）　　＋ｎ　＋
　　（Ｍ、　Ｎ）　］　＋ｒｐ　　［ｎ　（Ｍ、　Ｎ）
　＋ｎＩ　（Ｍ、Ｎ）］　ｌ　　　　　　　　　・・・
（１０）次に、電流方程式（３）、（４）を離散化する
。

この際、従来から知られているように、数値計算上の安
定性を確保するために、５ｅｈａｒｆ’ｅｔｔｅｒ−Ｇ
ｕｍｍｅｌの近似形を採用する。その結果、正孔、電子
の電流密度のＸ＋Ｖ成分は各々次のように与えられる。

Ｊ　、、（Ｍ’）−［ｑ　／　ｈ　、　（Ｍ’）］×　
［λｐａｌ　（Ｍ’）ｐ　（Ｍ　Ｎ）　＋λｐ−２（「
）ｐ　（Ｍ＋ｌ、Ｎ）　］・・・（１１−１） Ｊ　、、（Ｎ”）　−Ｃｑ　／　ｈ　ｙ　　（Ｎ　’　
）　］×［λ１１２１　（Ｎ’）ｐ　（Ｍ　Ｎ）十λ□
２　（Ｎ’）ｐ　（Ｍ、Ｎ＋１）　］（１１−２） Ｊ　、、（Ｍ’）＝　　［ｑ　／　ｈ　、　　（Ｍ’）
］×　［λｐ＋　　（Ｍ’）ｎ　（Ｍ　Ｎ）　　＋λ−
２（「）ｎ　（Ｍ＋ｌ、Ｎ）・・・　（１１−３） Ｊ　　−＝（Ｎ’）−［ｑ　　／　ｈ　　ｙ　　（Ｎ’
）コ×　［λｐ＋　　（Ｎ’）ｎ　（Ｍ　Ｎ）　　＋λ
ａｙ２（Ｎ’）ｎ　（Ｍ、Ｎ”ｌ）・・・　（１１−４
） ただし、上式のλは次式で与えられる。

λ　−−＋　　　ＣＭ’）−Ｃμ　ｐ　　　（Ｍ’）／
　　θ　（Ｍ　゛）コ　×［β（Ｍ　’）／　（１−ｅ
−β（Ｍ−１）　］λ−２ＣＭ’）−［μ、（Ｍ’）／
θ（Ｍ’）］Ｘ［β（Ｍ　’）／　（１−ｅ　′ｌ−’
　）　］］λ−７ｌ（Ｎ’）−［μ、（Ｎ’）、／θ（
Ｎ’）］ｘ［β（Ｎ　’）／　（１−ｅ　−”Ｎ−’）
　３λ−７２（Ｎ’）−［μ、（Ｎ’）／θ（Ｎ”）］
　Ｘ［β（Ｎ“）／　（１−ｅ””’　）コλ−＋　　
（Ｍ’）＝　［μ、（Ｍ’）／μ、（Ｍ’）］×λ−２
（Ｍ’） λ　。−２（Ｍ’）−Ｅ　μ　。　　（Ｍ’）／　μ　
、　　　（Ｍ’）コ×λｐ、Ｉ　　（Ｍ’） λｐ、Ｉ　　（Ｎ’）−［μ７（Ｎ　’）／　ｕ　ｐ　
　（Ｎ　’）］］ ］ ×λ−Ｆ２　　（Ｎ’） λ１□　（Ｎ’）−［μ、（Ｎ’）／μｐ　　（Ｎ’）
］×λｐ１　（Ｎ’） ・・・・・・　（１２） 以上により、原方程式（１）〜（５）の離散系への書き
替えが終了した。即ち、我々が解くべき方程式は、（８
−１）、（８−２）式に電流の式（１１−１）ないしく
１ｌ−４）、補助関係式（１２）と再結合の式（１０）
とを代入したものと、式（９）であり、これらは３個の
未知量ｐ１ｎ１ψを有する３個の方程式からなる。とこ
ろが、（８−１）、（８−２）式を見ると、これらは未
知量につき非線形な関係式を含むから１．このままでは
解を求めることが困難である。そこで、次に、これらの
式を未知量につきティラー展開し、２次以上の項を無視
することにより、未知量の変化分δｐ、δｎ、δψにつ
いての線形方程式を導出する。

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

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

Ｊ□（Ｍ′）叱Ｊ、、。（Ｍ’　） ＋　［ａＪ、、。（Ｍ’　）／ａｐ　（Ｍ、Ｎ）］×δ
ｐ　（Ｍ、Ｎ） ＋［ａ　Ｊｐ−”　　（Ｍ’　）　／ａ　ｐ　（Ｍ’＋
　１．　Ｎ）コ×δｐ　（Ｎ＋１．Ｎ） ＋　［ａＪ　１．’　　（Ｍ’　）／∂ψ（Ｍ、Ｎ）］
×δψ（Ｍ、Ｎ） ＋　［ｆｌｌｌＪ、、’　　（Ｍ’　）／’∂ψ（Ｎ＋
１．Ｎ）］×δψ（Ｎ＋１．Ｎ）　　　　　・・・（１
４−，１）Ｊ、、（Ｎ’）七Ｊ、、　　（Ｎ’） ＋　［ａ　Ｊｏｙ’　　（Ｎ’　）／ａ　ｐ（Ｍ、Ｎ）
コ×δｐ（Ｍ、Ｎ） ”　　［ａＪｐｙ’　　（Ｎ’　　）／９ｐ　　（Ｍ、
　　Ｎ＋１）］×　δ　ｐ　（Ｍ、　　Ｎ＋１） 十　　［９ＪＰＦ。　　（Ｎ’　　）　　／ａ　　ψ　
　（Ｍ、　　　Ｎ）　　］Ｘ　δ　ψ　（Ｍ、　　Ｎ） ＋　　［８Ｊ　、、’　　（Ｎ’　　）／ａ　　ψ　（
Ｍ、　　Ｎ＋１）］×δψ　（Ｍ、　　Ｎ＋１）　　　
　　　　　・・・　（１４−２）Ｊ、、（Ｍ’）　　た
Ｊ、、　　　（Ｍ’）＋　［９Ｊ−−（Ｍ’　　）／ａ
　　ｎ　　（Ｍ、　　Ｎ）］×　δ　ｎ　　（Ｍ、　　
Ｎ） ＋　［ａ　Ｊ、、’　　（Ｍ’　）／ａ　ｎ　（Ｎ＋１
．Ｎ）］×　δ　ｎ（Ｎ＋１．Ｎ） ＋　　［ａ　　Ｊ、、　　　（Ｍ′　＞　　／ａ　　ψ
　（Ｍ、　　Ｎ）］Ｘ　δ　ψ　（Ｍ、　　Ｎ） ＋　［ａ　Ｊ　、、’　　（Ｍ’　）／ａ　ψ（Ｍ＋　
１．Ｎ）コ×δψ（Ｎ＋１．Ｎ）　　　　　　・・・（
１４−３）Ｊ、、（Ｎ’）　　たＪ、、’　　（Ｎ’　
　）＋　　［１３Ｊ−ｙ　　　（Ｎ’　　）　　／９　
　ｎ　　（Ｍ、　　Ｎ）　　］×　δ　ｎ（Ｍ、Ｎ） ＋［ａ　Ｊ−ｙ’　　（Ｎ’　）／ａ　ｎ　（Ｍ、Ｎ＋
１）］Ｘδ　ｎ　　（Ｍ、　　Ｎ＋１） ＋［ａＪ−ｙ　　（Ｎ’・）／９ψ（Ｍ、Ｎ）コ×δψ
　（Ｍ、　　Ｎ） ＋　［８Ｊ、、’　　（Ｎ’　）／′∂ψ（Ｍ、Ｎ＋１
）］×δψ　（Ｍ、Ｎ＋１）　　　　　　　・・・　（
１４−４）ただし上つきゼロをもつ各員は未知量の試行
値ｐ　Ｏ、ｎ　Ｏ１ψ０に対応するものとする。

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

Ｕ（Ｍ、Ｎ）　−Ｕ’　（Ｍ、Ｎ）　十Ｕ、　’　（Ｍ
、Ｎ）　Ｘδｐ（Ｍ、Ｎ）　　＋　Ｕ　、　　’　　（
Ｍ、Ｎ）　　δｎ　（Ｍ、Ｎ）ただし Ｕ　、　　（Ｍ、Ｎ）　　−ａ　Ｕ　（Ｍ、Ｎ）　　／
　ａ　ｐ　（Ｍ、Ｎ）＝　　［ｎ（Ｍ、Ｎ）　　−ｒ　
　、、　Ｕ（間、Ｎ）］／（ｒ。

×［ｐ（Ｍ、Ｎ）　　＋　ｎ　＋　　（Ｍ、Ｎ）］＋　
ｒ　ｐ　　［ｎ　（Ｍ、Ｎ）＋　ｎ　＋　（Ｍ、Ｎ）１
１ Ｕ　、　（Ｍ、Ｎ）　＝　ｆ３　Ｕ　（Ｍ、Ｎ）　／　
ａ　ｎ　（Ｍ、Ｎ）−［ｐ　　（Ｍ、Ｎ）　　−τ　、
Ｕ（Ｍ、Ｎ）　　コ　／　　（τ　。

Ｘ　［ｐ　（Ｍ、Ｎ）　＋　ｎ　＋　（Ｍ、Ｎ）］＋　
ｒ　ｐ　　［ｎ　（Ｍ、Ｎ）＋　ｎ　＋　（Ｍ、Ｎ）］
）　　　　　　　　・・・（１５）ここで、（１３）、
（１４−１）、（１４−２）。

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

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

ａｚ−−［１／Ｑｈｙ　　（Ｎ）］ ｘａＪ□。（Ｎ　−１）／　ａ　ｐ（Ｍ、Ｎ−１）ａ１
２驕Ｏ ａ　＋ｉ−［１／　ｑ　ｈｙ　　　（Ｎ）　　］ｘ　ｅ
　　Ｊ　ｐｙ。（Ｎ　　−１）／　ｆａ　　ψ（Ｍ、Ｎ
−１）ａ２　＋−〇 ａ２□−−［１／　ｑ　ｈｙ　　　（Ｎ）　　］Ｘ　ａ
　　Ｊ　、。（Ｎ　−１）／　９　ｎ　（Ｍ、Ｎ−１）
ａ２ｉ−［１／ｑｈｙ　　　（Ｎ）］ Ｘ　ｅ　　Ｊ　、、ｏ　　（Ｎ　　−１）／　９　ψ（
Ｍ、Ｎ−１）ａ　３１−０ ａ３２″″Ｏ ａ　３ｉ−７＋　　（Ｍ、Ｎ） −１／［ｈ、　　　（Ｎ）ｈ、（Ｎ−１）］ｂ＋＋−［
１／ｈ−（Ｍ）］ Ｘ［ａＪ　□’　　（Ｍ　　）／　９　　ｐ　（Ｍ、Ｎ
）−９Ｊ　、、。（Ｍ　−１）／　ｆＪ　ｐ（Ｍ、Ｎ）
　　］＋［１／ｈ、　　　（Ｎ）］ Ｘ　　［９Ｊ　、、’　　（Ｎ　　）／ａ　　ｐ　（Ｍ
、Ｎ）−Ｂ　Ｊ　、、’　　（Ｎ　−ｆ）／　Ｅｌ　ｐ
　（Ｍ、Ｎ）　　］＋　ａ　Ｕ　’　　（Ｍ、Ｎ）　　
／　ａ　　ｐ　（Ｍ、Ｎ）ｂ　、２−　ｉｌｌ　Ｕ　’
　　（Ｍ、Ｎ）　　／　ｆ９　ｎ　（Ｍ、Ｎ）ｂＢ−［
１／ｑｈ、　　　（Ｍ）］ Ｘ　　［ａ　　Ｊ　、、ｏ　　（Ｍ　　）／９　　ψ（
Ｍ、Ｎ）−９Ｊ　、ｏ　　（Ｍ　　−１）／　ｉ３　　
ψ（Ｍ、Ｎ）　　］＋［１／ｑｈ、　　　（Ｍ）］ Ｘ　　Ｃｃ）　　Ｊ　、、’　　（Ｎ　’）／Ｅｌ　ψ
（Ｍ、Ｎ）−９Ｊ　、、’　　（Ｎ　　−１）／　ａ　
　ψ（Ｍ、Ｎ）　　３ｂ　、、−−ａ　Ｕ　０（Ｍ、Ｎ
）　　／　ａ　　ｐ（Ｍ、Ｎ）ｂ２□−［１／ｑｈ、　
　（Ｍ）コ Ｘ　　［ａ　　Ｊ　、、’　　（Ｍ　　）／　ａ　　ｎ
　（Ｍ、Ｎ）−ａ　Ｊ　−、’　（Ｍ　−１）／　ａ　
ｎ　（Ｍ、Ｎ）コ十［１／ｑｈ、　　（Ｎ）　］ ×［ａ　　Ｊ　−ｙ’　　（Ｍ　　）／　ａ　　ｎ　（
Ｍ、Ｎ）−ｆＪ　　Ｊ　−ｙ’　　（Ｎ　−１）／　ａ
　ｎ　（Ｍ、Ｎ）　　］−ａ　Ｕ　’　　（Ｍ、Ｎ）　
　／　ａ　　ｎ　（Ｍ、Ｎ）ｔｚｉ”　　［１／ｑ　ｈ
、　　　（Ｍ）］Ｘ　　［８Ｊ　、、ｏ　　（Ｎ　　）
／ａ　　ψ（Ｍ、Ｎ）−９Ｊ　、、’　　（Ｍ　−１）
／　ａ　　ψ（Ｍ、Ｎ）　　］＋［１／ｑｈ、　　（Ｎ
）コ ｘ　　［ａ　　Ｊ　、ｏ　　（Ｎ　　）／ａ　ψ（Ｍ、
Ｎ）−９Ｊ　、、ｏ　　（Ｎ　　−１）／　ａ　ψ（Ｍ
、Ｎ）　　］ｔ）ｉｔ−Ｑ／ε　　　　　　　ｌ）　、
　２ｍ　−ｑ　／εｂ３３”−［１／ｈ、　　　（Ｍ）
］　　［１／ｈ、（Ｍ　　）＋１／ｈ、（Ｍ−１）コ ［１／ｈ、　　（Ｎ）＋１／ｈ、（Ｎ−１）］］＋＋＋
−［１／　Ｑ　ｈ、　　　（Ｎ）　　］Ｘ　　［Ｆ３Ｊ
　ｐｙ’　　（Ｎ　）／　ａ　ｐ　（Ｍ、Ｎ＋１）　　
］］＋１□譚０ （＋ｉａ　　［１／ｑｈｙ　　　（Ｎ）］Ｘ　［ａ　Ｊ
　、Ｆ。（Ｎ　）／　９　ψ（Ｍ、Ｎ＋１）コＣ２，ｍ
Ｑ Ｃ２□−［１／ｑｈ、　　　（Ｎ）］ Ｘ　　［ａ　Ｊ　、、ｏ　　（Ｎ　）／　ａ　ｎ　（Ｍ
、Ｎ＋１）　　］Ｃｚｉｍ　　［１／ｑｈｙ　　　（Ｎ
）］Ｘ　　［ｆＪ　　Ｊ　、、ｏ　　（Ｎ　　）／９　
ψ（Ｍ、Ｎ＋１）　　］］＋３＋−０、Ｃ３２−０ ｅｉｉ−１／ｈｙ　　　（Ｎ）ｈ、（Ｎ　　）ｄ＋＋−
［１／ｑｈ−（Ｍ）３ Ｘ　　［ａ　　Ｊ　□ｏ　　（Ｍ　−１）／　ｉｌｌ　
　ｐ　（Ｍ−１，Ｎ）　　］ｄ　、２−　Ｏ ｄ＋ｉ＝　　　［１／ｑｈ、　　　（Ｍ）］ｘ　　［ａ
　Ｊ　、、ｏ　　（Ｍ　−１）／ａ　ψ（Ｍ−１，Ｎ）
ｄｚ＋−０ ｄ２２−−　　［１／ｑ　ｈ、　　　（Ｍ）］ｘ　［ａ
　Ｊ　、、’　（Ｍ　−１＞／　ａ　ｎ　（Ｍ−１，Ｎ
）ｄｚｉ−［１／ｑｈ、　　　（Ｍ）］ ｘ　［９Ｊ　、、’　（Ｍ　−１）／　ａ　ψ（Ｍ−１
，Ｎ）ｄｓ＋＝ｏ、ｄ３□−０ ｄｉｉ−１／ｈ−（Ｍ）ｈ、（Ｍ　−１）ｅ　ｚ＝　　
［１／　ｑ　ｈ−（Ｍ）　　］Ｘ　　［９Ｊ　□ｏ　　
（Ｍ　）／　ａ　ｐ（Ｎ＋１．Ｎ）　　］ｅ１□−Ｏ ］＋３−　　［１／ｑｈ、　　　（Ｍ）］Ｘ　　［’ａ
　Ｊ　ｅ、。（Ｍ　）／　９　ψ（Ｎ＋１．Ｎ）　　］
ｅ２１″″０ ｅｚｚ”［１／ｑｈ、　　　（Ｍ）］ Ｘ　　［ａ　Ｊ　＋ｔｔ’　　（Ｍ　）／　９　ｎ　（
Ｎ＋１．Ｎ）　　］ｅ２ｓ−［１／　ｑ　ｈ−（Ｍ）］ ｘ　　［ａ　Ｊ　ｓｗ’　　（Ｍ　）／ａ　ψ（Ｎ＋１
．Ｎ）　　］ｅ　３１　”　Ｏｒ　　ｅ　３２ｗ　０ｅ
３３−１／ｈ、　　　（Ｍ）ｈｘ　　（Ｍ　）ｆ　Ｉ−
Ｕ　’　（Ｍ、Ｎ）　十Ｕ　ｅ　’　（Ｍ、Ｎ）　ｐ’
　（Ｍ、Ｎ）十Ｕ、　’　（Ｍ、Ｎ）　ｎ　’　（Ｍ、
Ｎ）［１／ｑｈ、　　　（Ｍ）］ Ｘ　　［ａ　　Ｊ　＄　　’　　（Ｍ　−１）／ｌ　ψ
（Ｍ−１，Ｎ）　　］×　［ψ’　　（Ｍ−１，Ｎ）　
　−ψ’　　（Ｍ、Ｎ）　　］＋［１／ｑｈ、　　　（
Ｍ）］ Ｘ［ａＪ、’（Ｍ）／ｆ３　ψ（Ｍ、Ｎ）　　］×　［
ψ’　　（Ｍ、Ｎ）　　−ψ’　　（Ｎ＋１．Ｎ）　　
］［１／ｑｈ、　　　（Ｎ）］ Ｘ　　［９Ｊ　、　　’　　（Ｎ　　−１）／ａ　　ψ
（Ｍ、Ｎ−１）］Ｘ　［ψ’　　（Ｍ、Ｎ−１）　　−
ψ’　　（Ｍ、Ｎ）　　］＋　　Ｅｌ／ｑｈ、’　　（
Ｎ）］ ｘ　　［ａ　　Ｊ、　　’　　（Ｎ　　）／ａ　　ψ（
Ｍ、Ｎ）　　］×　［ψ’　　（Ｍ、Ｎ）　　−ψ’　
　（Ｍ、Ｎ＋１）］ｆ２−Ｕ’、（Ｍ、Ｎ）　　−Ｕｐ
　　’　　（Ｍ、Ｎ）　　ｐｏ（Ｍ、Ｎ）−Ｕ　＠　　
’　　（Ｍ、Ｎ）　　ｎ　’　　（Ｍ、Ｎ）［１／ｑｈ
、　　　（Ｍ）］ Ｘ　　［１３Ｊ　、　　’　　（Ｍ　−１）／ａ　ψ（
Ｍ−１，Ｎ）　　］×　［ψ’　　（Ｍ−１，Ｎ）　　
−ψ’　　（Ｍ、Ｎ）　　］＋［１／ｑｈ、　　　（Ｍ
）］ ｘ　−Ｃｆ３　Ｊ　、　’　（Ｍ　）／７３　ψ（Ｍ、
Ｎ）　］Ｘ［ψ’　（Ｍ、Ｎ）−ψ’　（Ｍ＋Ｉ、Ｎ）
コ［１／Ｑｈ、　　　（Ｎ）］ ｘ　　［９Ｊ　、　’　（Ｎ　−Ｄ、／　ｆ３　ψ（Ｍ
、Ｎ−１）×　［ψ’　（Ｍ、Ｎ−１）　　−ψ’　（
Ｍ、Ｎ）　］＋　　［１／ｑｈ、　　　（Ｎ）］ ｘ　［ａ　Ｊ　、　’　（Ｎ　）／ａ　１Ｍ、Ｎ）　］
×［ψ’　（Ｍ、Ｎ）−ψ’　（Ｍ、Ｎ＋１）コｆ３　
＝　　　（ｑ／ε）ｒ’　　（Ｍ、Ｎ）・・・　（１８
） 以上により、離散化した非線形方程式（８−１）。

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

式（１６）は、２次元配列された多数の分割点の１１分
に相当するものであるから、デバイス平面全体について
解を得るためには式（１６）を全分割点につき連立した
ものを解かねばならない。

いま仮にデバイス全体に対応する分割点番号が、１≦Ｍ
≦に、１≦ＮＩＬの範囲にあるものとすれば、式（１６
）を全体につき書き下した結果は第１０図に示すように
なる。但しここで簡単のため、ＡＴ　（Ｍ、Ｎ）をＡＭ
Ｎ、ｅｒ　　（Ｍ、Ｎ）をｅ　ＭＮ。

ＦＴ　　（Ｍ、Ｎ）をＦＭＮなどと書きかえた。

第１０図において％Ｂ１１など個々の行列は（３Ｘ　３
）のディメンションをもつから、全体の行列のディメン
ションは（３ＫＬ）Ｘ　（３ＫＬ）となる。いま標準的
なデバイスモデルとしてに−３０，Ｌ−４０を与えるも
のとすれば、全体の行列のディメンションは８．８０Ｑ
Ｘ　　３．８００−１２．９６０．０００−　１．２９
６Ｘ　１０　’の大きさとなる。

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

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

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

１≦Ｎ≦Ｌ］を改めて［ｐｏ（Ｍ、Ｎ）　、　　ｎ’　
（Ｍ、Ｎ）　。

ψ’　（Ｍ、Ｎ）　、　　１≦Ｍ≦に、１≦Ｎ≦Ｌ］と
置き換え、再度上述したのと同一の計算プロセスを実行
する。以下同様にして、必要回数だけ反復計算を実行す
れば、遂には非線形方程式（８−１）。

（８−２）、（９）を完全に満足する解が求まる。

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

この問題の解法のため、本発明は、従来法から離れて、
行列解法を適用することなしに所定の方程式の解を求め
、半導体デバイスの動作解析、所定の物理現象の解析を
行う方法及び装置を提供するものである。

［発明の構成］ （課題を解決するための手段） 本発明の第１の半導体デバイスの動作解析法は、半導体
デバイスモデリングにおける、電子、正孔の輸送方程式
およびポアソン方程式からなる連立方程式系をコンピュ
ータを用いて解析する方法において、前記連立方程式系
ｆ、−０，ｆ、−０゜ｆ、−０を、人工的時間微分項ｄ
ｐ／ｄｔ、ｄｎ／ｄｔ、ｄψ／ｄｔと感度係数λ２．λ
ｐ、λ。

とを含む下記（ａ　１）　、　　（ｂ　１）　、　　（
ｃ　１）式に置き換え、 ｄｐ／ｄｔ−一λｐｆｐ　　　　　　　”’（ａｌ）ｄ
　ｎ　／　ｄ　　ｔ　−λｐ、　ｆ、　　　　　　　　
−（ｂｌ）ｄψ／ｄｔ　−λｄｆ４　　　　　　　・・
・　（ｃ１）半導体デバイスの分割点配置（Ｍ、Ｎ）を
決定すると共に、感度係数λ２．λ０．λｐにも前記半
導体デバイスの物理的特性に合致した適正な空間位置依
存性を与えて、上記（ａｌ）、（ｂｌ）。

（ｃｌ）式を下記の（ｄ）、（ｅ）、（ｆ）式に変換し
、 ｄｐ（Ｍ、Ｎ）／ｄｔ−−λｐ　（Ｍ、Ｎ）　　ｆ　ｐ
（Ｍ、Ｎ）・・・（ｄ） ｄ　ｎ　（Ｍ、Ｎ）　／　ｄ　を−λｐ　（Ｍ、Ｎ）　
　ｆ　、　（Ｍ、Ｎ）・・・（ｅ） ｄψ（Ｍ、Ｎ）　／　ｄ　ｔ−λａ　（Ｍ、Ｎ）　ｆ　
ａ　（Ｍ、Ｎ）・・・（ｆ） ただし各式の右辺のｆ、、ｆ、、ｆ、は次式のとおり定
義し、 ｆ　、　（Ｍ、Ｎ）　−（１／ｑ）　　（［Ｊ　、、（
Ｍ　）−Ｊ、、（Ｍ−１）］／ｈ、　　（Ｍ）１＋　（
１／ｑ）（［Ｊ、、（Ｎ　） −Ｊ□（Ｎ　−１）］　　／ｈ、　　　（Ｍ）　　）＋
　Ｕ　（Ｍ、Ｎ）　　　　　・・・　（２２−１）ｆａ
　　（Ｍ、Ｎ）　　−（１／ｑ）　　Ｉ　　ＥＪ、、Ｃ
Ｍ　）−Ｊ　−（Ｍ　−１）］　　／　ｈ　、　　　（
Ｍ）＋　　（１／ｑ）　　（［Ｊ、、（Ｎ　　）−Ｊ□
（Ｎ　−１）］　　／ｈ、　　　（Ｍ）　　１−　Ｕ　
（Ｍ、Ｎ）　　　　・・・　（２２−２）ｆａ　　　１
Ｍ、Ｎ）　　　−［１／ｈ−（Ｍ）　　］×（［ψ（Ｎ
＋１．Ｎ）　　−φ（Ｍ、Ｎ）　　３　　／ｈ、　　（
Ｍ　）［ψ（Ｍ、Ｎ）　　−ψ（Ｍ−１，Ｎ）　　］　
　／　ｈ　、　　（Ｍ　−１））＋［１／ｈ、　　（Ｎ
）］ ×（［ψ（Ｍ、Ｎ＋１）　　−φ（Ｍ、Ｎ）　　］　　
／ｂ、　　（Ｎ　’）［ψ（Ｍ、Ｎ）　　−ψ（Ｍ、Ｎ
−１）　　］　　／　ｈ　、　　（Ｎ　−１）１＋　（
ｑ／ε）　　［ｒ’（Ｍ、Ｎ）　　＋ｐ（Ｍ、Ｎ）−ｎ
　（Ｍ、Ｎ）　　］ ・・・　（２２−３） 上記（ｄ）、（ｅ）、（ｆ）式を時間積分することによ
り、前記連立方程式系の解を求めることを特徴とする。

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

本発明の第２の半導体デバイスの動作解析法は、半導体
デバイスモデリングにおける、電位と空間電荷の関係を
表わすポアソン方程式をコンピュータにより解析する方
法において、前記ポアソン方程式を時間微分項ａφ／ｅ
ｔと感度係数λを含む下記（ａ２）式に置き換え、 ∂ψ／ａ　ｔ−Ａ　ｆｔ　　　　　　　−（ａ２）半導
体デバイスの分割点配置（Ｍ、Ｎ）を決定すると共に、
感度係数λにも前記半導体デバイスの物理的特性に合致
した適正な空間位置依存性を与えて、上記（ａ２）式を
下記の（ｂ２）式に゛変換し、 θψ（Ｍ、Ｎ）　　／　９　　を−λ（Ｍ、Ｎ）　　ｆ
　、　　０４．）ｌ）・・・　（ｂ　２） 上記（ｂ２）式を時間積分することにより、前記ポアソ
ン方程式の解を求めることを特徴とする。

本発明の第３の半導体デバイスの動作“解析法は、最小
自乗法に基く半導体デバイスの等倍回路定数決定問題を
コンピュータにより解析する方法において、前記等価回
路定数決定問題を時間微分項ｄｑ／ｄｔと感度係数λと
を含む下記（ａ３）式％式％ （３） 感度係数λに個々の等価回路定数の有する物理的特性に
適した値を与えて、前記（ａ３）式を時間積分すること
により、前記等価回路定数決定問題の解を求めることを
特徴とする。

本発明の第４の所定の物理現象の解析法は、所定の物理
現象を表わす単一の方程式または連立方程式系ｆ　（ｘ
）−０をコンピュータにより解析するために、前記方程
式または方程式系を時間微分項ｄ　ｘ　／　ｄ　ｔと感
度係数λとを含む下記（ａ４）式に置き換え、 ｄ　ｘ　／　ｄ　を＝λｆ（ｘ）　　　　　　−＝Ｃａ
４）感度係数λに前記物理現象に適した値を与えて、前
記（ａ４）式を時間積分することにより、前記行列方程
式の解を求めることを特徴とする。

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

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

本発明の第１の解析装置は、半導体デバイスモデリング
における、電子、正孔の輸送方程式およびポアソン方程
式からなる連立方程式系を時間微分項と感度係数とを含
む微分方程式ｄｐ、”ｃｉｔ−一λｐｆ、、ｄｎ／ｄｔ
−λ。ｆ、、ｄψ／ｄｔ−λａｃｍに置き換えて解析す
る装置において、前記装置は、単数または複数のｍ個の
論理セルと、前記各論理セルに結合し調整可能な増幅利
得λ。

を有するｍ個のアンプまたはこれと等価な乗算器と、前
記各論理セルに対してネットワーク状に設けられ前記装
置内の過度現象が定常状態になったことを検知する検出
手段とからなることを特徴とする。

本発明の第２の解析装置は、半導体デバイスモデリング
における電位と空間電荷の関係を表わすポアソン方程式
を微分方程式ｄψ／ｄｔ＝λｆ。

に置き換えて解析する装置において、前記装置は、単数
または複数のｍ個の論理セルと、前記各論理セルに結合
し調整可能な増幅利得λｐを有するｍ個のアンプまたは
これと等価な乗算器と、前記各論理セルに対してネット
ワーク状に設けられ前記装置内の過度現象が定常状態に
なったことを検知する検出手段とからなることを特徴と
する。

本発明の第３の解析装置は、最小自乗法に基く半導体デ
バイスの等倍回路定数決定問題を微分方程式ｄｑ／ｄｔ
−一λ（Ａｑ−ｂ）に置き換えて解析する装置において
、前記装置は、単数または複数のｍ個の論理セルと、前
記各論理セルに結合し調整可能な増幅利得λｐを有する
ｍ個のアンプまたはこれと等価な乗算器と、前記各論理
セルに対してネットワーク状に設けられ前記装置内の過
度現象が定常状態になったことを検知する検出手段とか
らなることを特徴とする。

本発明の第４の解析装置は、所定の物理現象を表わす単
一の方程式または連立方程式系ｆ　（ｘ）−〇の解を求
める計算を微分方程式ｄ　Ｘ／ｄ　ｔ　＝λｆ　（ｘ）
の積分問題に置き換えて実行するように組み立てられた
解析装置が、時間軸上の積分プロセスを実行するアナロ
グ計算実行部と、計算途上で係数行列λの大きさを個別
に調節、適正化する手段とを有することを特徴とする。

（作　用） 本発明の方法では、半導体デバイスモデリングにおける
電子、正孔の輸送方程式およびポアソン方程式からなる
連立方程式系の解法、半導体デバイスモデリングにおけ
る電位と空間電荷の関係を表わすポアソン方程式の解法
、最小自乗法に基く半導体デバイスの等価回路定数決定
問題の解法等において、上記各式を、時間微分項と感度
係数λとを有する微分方程式に置き換え、感度係数λの
値を解くべき式の物理的特性に合わせて設定するので、
対象となる方程式がさらに大形化しても、行列解法を高
速に実行できる。

本発明の装置は、式（ａ４）の独立変数ｔを時間とみな
し、未知ｆｆ１ｘは式（ａ４）を満足する電気回路網の
所定端子の電圧で表わされるものとする。この回路網の
過渡現象が時間と共に進行して定常状態が実現すれば、
それがすなわち式（ａ４）の解を与えるように装置全体
が構成される。

一般に、従来法に基づくディジタル計算では、行列方程
式Ａｘ＝ｂの解の計算において、消去法を実行する際は
、形式上ｘ−Ａ−’ｂと書かれる、所謂逆行列計算ない
しは実質上これと等価な計算が含まれ、これを実行する
ために多大の計算時間が消費される。

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

順方向の計算が有利なことは、特に行列Ａが多数のゼロ
要素を含む場合において著るしい。この場合かけ算Ａｘ
’の実行に際してＡの非ゼロ要素のみを考慮すればよい
からである。本計算法が行列方程式Ａｘ−ｂの真の解を
与えるためには、ひとつの条件がある。それは、微分方
程式（ａ４）の、時間軸上の積分計算が十分長く実行さ
れた場合、未知量Ｘが定常値に到達することである。こ
れを式で書けば、 分な時間が経過したあとは、微分方程式（ａ４）の左辺
はゼロとなる。従って右辺もゼロでなければならないこ
とから、求まったＸはＡｘｍｂを満足する。

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

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

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

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

次にλの適正化法を論じる。λは前述したとおり対角行
列であり、ｎ個の要素λ０．λ２．・・・・・・λイを
もつ。すなわち λ−旧ａｇ　（λ１．λ２．・・・λ。）・・・（１つ
） である。ここで式（２〉において、未知ｆｌＸ＋、ｊ−
１，２，・・・、ｎの各々につき、装置系全体が安定動
作し得るような最大の時間微分量絶対値ｄ　Ｘ　Ｉ／　
ｄ　ｔ　ｌ　ｓａｗが存在するはずであるから、λの成
分λｐは、るこれからきまる限界値を超えてはならない
。これを式で表わせば、 Ｘｌ　１≦（１／ＩΣａ＋ＩＸ１　　ｂｔ　　Ｉ）Ｘｌ
　ｄ　ｘ　＋／　ｄ　ｔ　Ｉ　ｍａｗ　。

ｉ−１，２，・・・、ｎ ・・・（２０） となる。式（２０）の右辺の大きさがどの程度になるか
、は、問題の性質により異なるので、一般論の範囲では
これ以上具体的に論じることは困難と思われる。行列状
数学によればこの種の判定条件は行列の固有値と関連づ
けられるのが常であるが、実際問題を解くことと行列の
固有値を求めることとはほぼ同質かつ同程度の困難さを
もつ問題であるため、実用上有効な判定条件を与えるこ
とは希である。そこで、実際に本発明の装置を組立てる
に際しては、後述のようにｎ個の論理セルの各々に付随
したアンプを設け、その各々が増幅利得λ１をもつよう
構成する。各アンプは個々の実際問題に十分対応可能な
程度の十分大きな増幅利得をもつように設計し、利得調
整用のつまみまたはコントローラを具えるものとする。

（実施例） 以下、図面を参照して本発明を説明する。

第１図は、本発明の第１実施例に係る半導体デバイスの
動作解析法を説明するためのフローチャートである。

本発明では、方程式（８−１）、（８−２）。

（９）の解を求めるにあたり、これらを時間微分項を含
む、つぎの形におきかえる。

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

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

ｐ　’　（Ｍ、Ｎ）　−ｐ　０（Ｍ、Ｎ）−δｔλｐ　
（Ｍ、Ｎ）　　ｆ　、　。（Ｍ、Ｎ）・・・（２３−１
） ｎ　’　　（Ｍ、Ｎ）　　−ｎ　’　　（Ｍ、Ｎ）＋δ
　ｔ　λ　ｎ　（Ｍ、Ｎ）　　ｆ　、　　　（Ｍ、Ｎ）
・・・　（２３−２） ψ’　　（Ｍ、Ｎ）　　−ψ’　　（Ｍ、Ｎ）＋δ　ｔ
　λａ　　（Ｍ、Ｎ）　　ｆ　ａ　　’　　（Ｍ、Ｎ）
・・・　（２３−３） ただし、１≦Ｍ≦に、１≦Ｎ≦Ｌとする。

以下、第１図のフローチャートを参照して説明する。本
方法においては、まず計算をスタートするにあたって、
ステップａ１において、半導体デバイスの構造、不純物
データをコンピュータのメモリに入力する。次にステッ
プｂ１において、半導体デバイスの分割点配置を決定す
ると共に、ステップＣ１でバイアス電圧Ｖを人力する。

次に基本変量ｐ（Ｍ、Ｎ）、ｎ　（Ｍ、Ｎ）、　ψ（Ｍ
、Ｎ）の初期試行値１）’　　（Ｍ、Ｎ）、ｎ’　　（
Ｍ、Ｎ）。

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

次にステップｆ１において、式（２３−１）〜（２３−
３）に従って各式の右辺の値をすべての点１≦Ｍ≦に、
１≦Ｎ≦Ｌについて求め、これらをｐ’　　（Ｍ、Ｎ）
、ｎ’　　（Ｍ、Ｎ）。

ψｌ　　（Ｍ、Ｎ）；１≦Ｍ≦に、１≦Ｎ≦Ｌに等しい
とおく。これによって時間軸上の積分が１ステップ進行
したことになる。

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

もしデバイス空間の中に１点でも誤差限界より大きな変
化を生じた点があれば、ステップｈ１において、先に求
めた修正値ｐｌ　　（Ｍ、Ｎ）。

ｎ’　　ＣＭ、Ｎ）、　　ψｌ　　（Ｍ、Ｎ）を新な試
行値とみなし、これらをＩ”　　ＣＭ　＋　　Ｎ　）　
、ｎ　０（Ｍ　。

Ｎ）、ψ０　（Ｍ、Ｎ）に等しいとおき、更に時刻１１
をｔｏとして、前述と同じ要領で更に時間軸上の積分計
算を行う。

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

（８−２）および（９）の解が求まったことに他ならな
い。

以上が本実施例の概略である。従来法とくらべた場合の
特徴は、■本実施例には非線形方程式の線形化プロセス
が含まれないこと、■従来法の説明の項で示した式（１
６）及びこれから派生した第１０図のような大形サイズ
の行列・ベクトル方程式の解の計算を必要とせず、従っ
て大形メモリ容量を用いた複雑な計算を実行することが
なく、単に式（２３−１）〜（２３−３）により各式の
右辺の値を求めて、基本変量の修正値を決めるだけでよ
いこと、■調整因子としての係数λｐ　（Ｍ、Ｎ）　、
　　λ。（Ｍ、Ｎ）およびλ＃　（Ｍ、Ｎ）　−一以下
これらを単に感度係数と呼ぶｍ−が、デバイス空間の各
点につき異なった値をとり得るため、これらに物理的特
性に合致した適正値を与えることにより、解の収束が効
率よく実現できること、■本実施例の式（２１−１）〜
（２１−３）の解を計算するプロセスは時間軸上の数値
積分計算に他ならない。この計算を通常の如く大形計算
機を用いたディジタル計算により実行することは勿論可
能であるが、他に、電子回路を応用することにより、そ
の回路の時間応答（過渡現象）が定常状態に達したなら
ば、その状態がすなわち所定の方程式の解に対応するよ
うな装置を構築することが可能である。この点について
は以下に本発明の装置として説明する。

以上、第１実施例の方法は、行列解法に基づ〈従来法と
は全く異なり、デバイスモデリングに新な分野を開拓す
る可能性をもつものである。

次に第２図を参照して、本発明の第２実施例に係る半導
体デバイスの動作解析法を説明する。

デバイスモデリングにおいて、たとえばｐｎダイオード
のｐ側をゼロ電位、ｎ側を正電位Ｖに保つ場合、ｐｎ接
合の界面近傍に空乏層が形成される。この際、電位Ｖが
ダイオードの耐圧以下であれば、ダイオードに流れる電
流はほとんど無視できる程度に小さい。このような場合
、基本方程式（１）〜（５）のうちボアリン方程式（５
）だけを解けばデバイス全体の電位ψの分布が正確に求
まることが知られている（例えば、Ｋｕｒａｔａ。

’ＮｕＩＩｅｒｉｃａｌ　Ａｎａｌｙｓｉｓ　ｆ’ｏｒ
　Ｓｅｍ１ｃｏｎｄｕｃｔｏｒＤｅｖ１ｃｅｓ、　Ｈｅ
ａｔｈ　、　１９８２．第８章、又は、倉出、「バイポ
ーラトランジスタの動作理論」、近代科学社、昭和５５
年３．２節および５．１節参照）。この際正孔密度およ
び電子密度ｐ、ｎは各々次式で与えられる。

ｐ＝ｎｔ　　ｅ−”　　　　　　　　　　（２４１）ｎ
　＝ｎ　ｒ　　ｅ”’−”　　　　　　　　（２４２）
これらを式（５）に代入すれば、２次元空間の場合次の
方程式が得られる。

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

これに対し本発明の方法では、時間微分項と感度係数λ
を含むっぎの形の方程式から出発する。

ａ　　ψ（Ｍ、Ｎ）　　／　９　　ｔ　　−Ａ　（Ｍ、
Ｎ）ｘ　　（［１／ｈ、　　’　　（Ｍ）　　］Ｘ（［
ψ　（Ｎ＋１．Ｎ）　　−ψ　（Ｍ、Ｎ）　　コ　／ｈ
、（Ｍ’）［ψ（Ｍ、　Ｎ）　　−ψ（Ｍ−１，Ｎ）　
　］　　／　ｈ　、　　（Ｍ’　　−１）　　１十［ｉ
　／ｈ　ｙ　’　（Ｎ）　］ Ｘ（［ψ（Ｍ、Ｎ＋１）　　−ψ（Ｍ、Ｎ）　　］　　
／　ｈ　、　　（Ｎ’　　）［ψ（Ｍ、Ｎ）　　−ψ（
Ｍ、Ｎ−１）　　］　　／　ｈ　、　　（Ｍ’　　−１
）　　１＋ρ（Ｍ、Ｎ）／ε）　　　　　　　　　　　
　・・・　（２８）この方程式を時間軸上に沿って積分
する手順は、すでに第１実施例の説明で述べたのと同一
要領で実行される。そのためのフローチャートを第２図
に示す。即ち、本方法においては、まず計算をスタート
するにあたって、ステップａ２において、半導体デバイ
スの構造、不純物データをコンピュータのメモリに入力
する。次にステップｂ２において、半導体デバイスの分
割点配置を決定すると共に、ステップｃ２でバイアス電
圧Ｖを人力する。

次に電位の初期試行値ψ０をステップｄ２で与え、ステ
ップｅ２において時間積分開始のための初期時間を日を
定める。

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

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

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

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

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

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

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

この反面、もし上記の変化の絶対値が大きすぎると、計
算過程が不安定となり、時間ｔの経過に伴ってψ（Ｍ、
Ｎ）の振動または発散を生じるから、真の解が求まらな
い。

感度係数λ（Ｍ、Ｎ）は、このような変量変化を適度に
生じさせるための調整因子である。特に、本例のように
半導体デバイスのポアッソン方程式で、かつ式（２９）
に示した２次元空間モデルの場合、上記の条件を満たす
ようなλ（Ｍ、Ｎ）として、次の式を用いればよいこと
が知られている。即ち、まず式（２９〉に関連して、次
の諸量を定義する。

Ｈ−１／　　［ｈ　　、　　（１４’−Ｄｈ　　、　　
°（Ｍ）コ＋　１／　　［ｈ、　　（Ｍ’）ｈ、　　°
００コ＋　１　／　［ｈ　、　（Ｎ’−１）ｈ　、　°
（Ｎ）］＋　１　／　［ｈ　、　（Ｎ’）ｈ　、　’（
Ｎ）］　　　・・・（３０−１）α海（ｑθｎ＋／ε）
× ［ｅｘｐ（−θψ’　（Ｍ、Ｎ）　１ ＋ｅｘｐ　　（θ（ψ０（Ｍ、Ｎ）　−Ｖ）　ｌ　］・
・・（３０−２） 以上により λ（Ｍ、Ｎ）　　・δｔ≦２／　（２Ｈ＋α）　・・・
（３０−３）とする。

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

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

次に第３図を参照して、第２実施例の方法を実行する装
置について説明する。このような装置を組み立てるには
、Ｎ７個の論理セルの各々に付随したアンプが増幅利得
λ１を有するように構成する。ここでＮＴは、連続系に
属する原方程式をコンピュータを用いて解くために、差
分法、有限要素法などの近似解法に基づき、原方程式を
離散系に置き換えた場合のメツシュの分割点の総数であ
る。各アンプは個々の実際問題に十分対応可能な程度の
大きな増幅利得を有するように設計し、利得調整用のつ
まみ又はコントローラを備えるものとする。ただし、問
題により増幅利得は１より小さいこともあり得るので、
この際には、アンプは実際上減衰器となるが、このよう
な場合も包括されるものとする。或いはまた、アンプと
等価な装置として、入力量と増幅利得の積を求めるため
の乗算器を用いてもよく、この場合も包括されるものと
する。

参照符号１〜１２は論理セルであり、非線形な入出力応
答を有する。この実施例の場合には、ＮＴは１２である
。各セルの入力端子にはアンプａｌ、ａ２．・・・ａ１
２がそれぞれ設けられている。

また、各アンプと各セルとのノードには、抵抗と容量と
からなる共振回路ｂｌ、ｂ２．・・・ｂｌ２が設けられ
ている。各アンプはλに相当する可変利得を有する。入
力電流源ｃｌ、ｃ２．・・・ｃ１２からは電流が各セル
に対して人力され、この入力端子が、本発明の装置内に
おいて定常状態になったことが、各セルに対してネトワ
ーク状に設けられた電流検出手段１１．２１．・・・１
２１２により検出される。

この検出結果はアドレスユニット１００によりアドレス
指定され、Ｄ／Ａコンバータ１０１でディジタル信号に
変換されて回路１０３に送られる。

回路１０３内には微分器、乗算器、′加算器が内蔵され
ていて、所定の演算を実行し、測定データを回路１０４
に出力する。

以上の実施例では、メツシュ分割点の総数Ｎ７に等しい
Ｎ７個のセルを用いた例を示したが、スィーパを使用す
ると、ＮＴ以下のｍ個のセルを用いて本発明の解析装置
を構成でき、上記実施例と同様な動作をさせることがで
きる。

次に第４Ａ、４Ｂ、５図を参照して、本発明の第３実施
例に係る半導体デバイスの等価回路定数の最小自乗法に
よる決定問題を説明する。いま、４個のｙパラメータ（
これらは全て複素アドミッタンスである）　ｙｌｌ＋　
　ｙ　Ｉ□＋　　Ｖ２＋＋　　Ｙ２２が、Ｎ個の周波数
ｆｋ　　（ｋ−１，２，、−、、Ｎ）１．：おイテＡ［
ｌＩ定されており、これをＭ個の等価回路定数ｑ＋　＋
４２１　・・・＋　　Ｑｕを含む等価回路（第４Ａ、４
Ｂ図参照）によりｆｉｔｔｌｎｇを試みる場合、つぎの
行列方程式が導かれる。

Ａｑ−ｂ、Ａ−（ａ＋１）、ｂ−（ｂｌ）・・・（３１
−１） ａ＋１−　Σ　　［α　　　ｊｓｋｌ　α　ｊ＋＋＋ｋ
ｌ＋　α　ｊｌ、ｌｋｌ　α　　　Ｉｍｋｌコ／１η７
ｍｋＢ （３１−２） ｍ−Σ　　［β　　７＋ｍｋ　　ａｌｍｋ（３１ ３） ただし上記において （Ｘｉｍｚ−Ｆ）ｙｌｍｋ　／９　ｑ＋　　　＋＋＋　
　（３２１）βｊ　＋ｗｋ　−Ｙ　ｐ　、、ｋ　　−η
ｊｍｋ−Σ　α、□ｑ　、　Ｏ・・・（３２−２）ここ
で、＊印は共役複素数を表わし、Ｎ＋　ｍは１゜２のい
ずれかの値を取るものとする。

上記行列方程式を、本発明の解法に基づきｄｑ／ｄｔ−
一λ（Ａｑ−ｂ）のように置き換える。この方程式は第
５図のフローチャートに基づき解かれる。ステップａ３
において、ｙノくラメータのＡｌ１１定データニηｌｌ
ｋ　＋　　η１２ｋ　＋　　η２１に一η２２ｉ　　；
ｋ　−１＋　　２　＋　・・・Ｎを読取る。次にステッ
プｂ３において未知量の初期試行値ｑ１Ｑ２　１　・・
・Ｉ　　ｑＭ’を与える。ステップＣ３において初期時
間と時間間隔を与え時間積分を開始する。ステップｄ３
において、ｙパラメータの理論値ｙ、、、とｙパラメー
タの測定値η、−とを用いてα、□１．β、１を求め、
ステップｅ３でａｌｌ＋ｂ１を求める。

次にステップｆ３において、ｑ＋　’をｉ−１゜２、・
・・１Ｍに対して計算する。ｑ　、　ｌは次式で与えら
れる。

ｑ＋’−ｑ＋　　　−６１−λ＋　　（ＡＱＯ−ｂ）・
・（３３） これによって時間軸上の積分が１ステ・ンプ進行したこ
とになる。

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

もし絶対値１ｑ、ｌ　　ｑ＋　　　ｌが１つでも誤差限
界より大きいときには、ステップｈ３において、先に求
めた修正値ｑ　、　ｌを新な試行値とみなし、これらを
ｑ、。に等しいとおき、更に時刻ｔＩをｔｏとして、前
述と同じ要領で更に時間軸上の積分計算を行う。

もしすべてのｉで変化量の絶対値が誤差限界を下まわっ
たならば、Ｍ個の等価回路定数ｑ１Ｑ２＋　・・・＋　
　ＱＭが得られたものとみなして５１算を終了する。そ
の結果得られた解は、ｄ　ｑ／ｄ　ｔ　−−λ（Ａｑ−
ｂ）において、式の左辺の時間微分項がゼロになるよう
なものであるから、方程式Ａｑ−ｂの解が求まったこと
に他ならない。

次に、トランジスタの等価回路定数決定問題につき、具
体的に実施した数値計算の結果を説明する。まず図４Ａ
に示したバイポーラトランジスタのエミッタ接地等価回
路を、同図４Ｂのごとく３個の４端子網の合成とみなせ
ば、ｙパラメータの計算を容易に行うことができる。途
中経過を省略して第４Ａ図、第４Ｂ図の真性部分１００
０のｙパラメータｙ、の計算結果を示せば、っぎのよう
になる。

Ｙ　＋＋＋　−（ｇ　ｂ　　［ｇ　＠　ｇ　ｍｅ　（１
−α）−ω２　ＣａＣｃｌコ ＋ｊωｇｂ　　（ｃｏ　ｇｓｅ＋ｃｃ１ｇ＊　＋Ｃｃ＋
ｇｓ、）／Δｙ１　　　　　　　　　　・・・（３４−
１）ｙ１□、　　ｗｍ、　　　［Ｊ　　（ＬＩ　Ｃｃｌ
ｇ　’×（ｇ、十ｇ、＋ｊ（Ｊ）ｃ、）］／Δｙｌ・・
・　（３４−２） Ｙ２１１　−　　Ｃｇｂ　　（αｇ、ｇ、、＋ω２　Ｃ
ａＣｃｌ）−ｊωＣｇ＋ｇ＋　　（ｇ、＋ｇｏ、）］　
　／ΔｙＩ・・・　（３４−３） ｙ２２１　　＝　　［−ω２Ｃａ　　Ｃｃｌ　（ｇｏ　
　＋　ｇ　＊Ｊ＋ｊωＣａｒ　（ｇｂ　　ｇｏ　　”ｇ
ｏ　　ｇ＊ｍ＋　ｇｃｅｇｂ　　）］／Δｙ−・・・　
（３４−４） ΔＹＩ−ｇｂｇ・　＋ｇｂｇ− ＋ｇ−ｇ＠、（１−α）−ω２　ｃ、　　ｃ。

十ｊω［Ｃ，（ｇ・・＋ｇｂ） 十Ｃ，Ｉ（ｇ、十ｇ０．）コ ・・・　（３４−５） 続いて４端子網の合成則により、真性部分１０００と外
部コレクタ容Ｍ１ｏｏ１を合成したもののｙパラメータ
ｙ１′は次式で与えられる。

ｙ＋＋−Ｖ＋＋＋＋ｊωＣ１゜・・・（３５−１）ＶＩ
２”−）’＋２１ｊωＣ６゜・・・（３５−２）’Ｗｚ
＋’　−Ｙ２ｚ　−ｊωｃｃ、−＝（３５−３）Ｙ　２
２”　＝　Ｖ　２２１　＋　ｊωＣｃ、−・・（３５−
４）更に合成則を用いてコレクタコンダクタンスｇ、を
含む全体すなわち１０００＆１００１＆１００２のｙパ
ラメータを計算する。途中を省略して結果だけ示せば、
次のようになる。

ｙｚ−（ｇゎＹｚ　　＋Δｙ”） ／　（ｇｇ　＋’）’２２”　）　　・・・（３６−１
）Ｙ　１２−　（ｇｃ　Ｙ　１２”　）　／　（ｇ・＋
Ｖ２２”）・・・（３６−２） Ｙ２＋−Ｃｇｔ　）’２１”　）　／　（ｇｅ　＋ｙ２
２’　）・・・（３６−３） Ｖ２２−　（ｇｃ　Ｙ２２”　）　／　（ｇｃ　＋Ｙ２
ｚ”　）・・・（３６−４） Δｙ”　−ｙ目　ｙ２２　３／１２　　Ｙ２１・・・（
３６−５） 第４Ａ図、第４Ｂ図において決定すべき回路定数は、ｇ
　ＩＩ　ｌ　　ｇ　＃＃＋　　ｇ　ｂ　ｌ　　ｇ　ｔ　
ｌ　　α、　　ｃ、、　　ｃｃｌ。

ｃ６゜の８個である。Ｎｔｔｉｎｇの指標となる目的関
数は、 Ｅ−Σ　ｌｙｎ□−η、□　１２／１η、□　１２Ｊ、
＋＋＋、　　叡 ・・・（３７） で与えられる。これはｙパラメータの理論値ｙ、＝１と
測定値ηＩ□のｒｆｔｌｌｌｌｇの相対誤差を表わす。

なおｌ、ｍは１，２のいずれがをとり、ｋは測定周波数
の番号に−１，２，・・・、Ｎを表わす。

次に計算実行上の便利のため、８個の未知量ｐ１〜ｐ８
および各々の規格化単位ｒ１〜ｒ８を下記のとおり定義
する。なお規格化の定義により、ｑ＋　−ｒ＋　ｐ＋　
、１−１．２．−８なることに注意する。

ｐ　、　ｍ　ｇ・　（１−α）ｌ　「１膳１ｐ２−αｇ
＊＋　　　　　　ｒ２■１ ｐｓ　−ｇ＊＊＋　　　　　　　　ｒ３ｍｌｐ４″ｗｇ
ｂ、　　　　　　　ｒ４−１ｐｓ　−ｇｃ　＊　　　　
　　　　ｍｓ　−１ｐ６−Ｃ・ｒ　　　　　　　　ｒ６
−１０−１２ｐｔ　”　Ｃｃｌ＋　　　　　　　　ｒ　
７■１０−＋２ｐ８　”　Ｃｃ　＋＋＋　　　　　　　
ｊ　ｓ　、ｌ　Ｑ　−１２ｙパラメータのｆｉｔｔｉｎ
ｇを行うには、その目標となる測定データが必要である
。本来は実際のデバイスの生データを用いて計算を行う
べきであるが、性質のよく分ったデータを用いたいとい
う要求から、回路定数ｐ１〜ｐ８に対し適当な値（以下
これを真値と呼ぶ）を与えて生成したｙバラメ−夕をも
って、人工的な「測定データ」とした。

以下の検討に用いた回路定数の真値はつぎの通りである
。

ｇｏ　−０，６ｓｉｅＩｌｅｎｓ　　　ｐ＋　−ｇｏ　
　（１−ｆｆ）−０，０１２ｓｌｅｉｅｎｓ α謬０．９８　　　　　　ｆ）２−ｇｏ　α＝　０．５
８８ｓｉｅａ＋ｃｎｓ ｇ　、、−０，０８ｓｉｅｍｃｎｓ　　−ｐ　　ｖｇ　
ｂ−０，０４ｓｌｅＩｌｌｅｎｓ　−ｐａｇ　ｔ　　−
０，１５ｓｌｅｍｅｎｓ　　−ｐ　ｓｃ　　、　　＝２
０Ｏｆ’Ｐ　　　　　　　−ｐｂＣｏ−３１Ｆ　　　　
−Ｉ−ｐ７ ｃ　　ｃａ−ｍ５０ｆ’Ｐ　　　　　　　　−ｐｇ「測
定」周波数は、０．１ＧＨｚを起点として０．５ＧＨｚ
毎に４０点を与える。その結果、最終値は１９．６ＧＩ
ｌｚとなる。

計算をスタートするため、回路定数の初期試行値として
、つぎの数値を与える。

ｇ　、　−０，７ｓｌｅＩＩｅｎｓ　　　ｐ　＋　−０
，０１４ｓ１ｅｍｅｎｓａ　−０，９８ｐ　２−０．６
８８ｓｌｅｉｅｎｓｇ　、、ｍ　Ｏ，１ｓｉｅａ＋ｅｎ
ｓ　　　−ｐ　３ｇ　ｈ　　−０，２ＳＩｅｆｆｉｅｎ
Ｓ　　　−９４ｇ　ｃ　　”　０．２ｓｌｅｍｅｎｓ　
　　−ｐ　９ｃ　ａ　−３０Ｏｆ’Ｐ　　　　−ｐ　６
Ｃ、＋ｍ　５　ｆＰ　　　　　−９７ Ｃｃ　ｏ　＝　８Ｏｆ’Ｆ　　　　　　　−ｐ　ｓ以上
に与えた真値と試行値に対応する４個のｙパラメータを
複素平面上にプロットして示すと、第６Ａ図〜第６Ｄ図
のようになる。測定データは、前述の通り人工的に生成
したものだから、生データに特有のノイジーなばらつき
を含まないものの、−見して分かるように同一周波数で
の真値と試行値のずれはかなり大きい。ずれの程度を定
量化するため、目的関数（３７）の平方根をとれば、ｙ
パラメータ全体の相対誤差ベクトルの大きさが分かる。

すなわち Ｒ，Ｅ、　＝　［、）％、、ｌ　ＶＪ＋ｍ＊　　７７＊
＊ｍ　ｌ　’／１η、□ｌ　２　、Ｉ　ｌ／２ ・・・（３８） とする。第６Ａ図〜第６Ｄ図の例は、Ｒ，Ｅ、−３，１
９すなわち３１９％である。望ましいＨｉｔｉｎｇの条
件を、大まかに誤差１％以下とすれば、上記の初期誤差
は非常に大きい。

時間軸上の積分を行うための時間きざみは、律に１０−
’５−Ｉｎｓとする。

実際に行なった、計算機による５回のテストランをまと
めて第７図に示す。最初のジョブＳ”は、一番左の欄に
示したλｐ〜λＢを用いて時間ステップ数−５００回の
計算を行った結果、相対誤差（Ｒ，Ｅ、　）−４，１２
％となる。これは初期誤差３１９％にくらべれば大幅な
改善であるが、まだ十分よいｆ’ｌｔｔｌｎｇとはいい
難い。そこで“Ｓ”の途中経過を参考にしながら、ｐｂ
　　（−ｃａ）の変化を大きく、ｐ７の変化を小さくし
た“Ｗｏを実行すると、Ｒ２Ｈ，−２，４７％と改善す
る。つぎの“Ｉ”では、“Ｗ”の最終結果を初期条件と
して同一のλで５００ステップ進行すると、さらに改善
してＲ，Ｅ、−１，０７％となる。以下同様にして“Ｑ
”では、“Ｗ”と同じ初期条件からλ７−５、Ｅ３　（
旧位−５，Ｅ２）でランすると、Ｒ，Ｅ、−０，２１０
％というすぐれた結果が得られる。このように−殻内傾
向として、まず第一回のランで誤差１桁％程度のほどよ
い結果を得たあと、その最終結果からスタートして、前
回と同じか、または大きめのλを用いて第二回目のラン
を実行すると、さらによい結果が得られる。

最後の′Ｕ”は“Ｑ”の最終値からスタートしたもので
、時間ステップ１９２回で最小誤差０．１％に到達する
。これが一連の計算で得た最善の結果であり、得られた
ｙパラメータを第６Ａ図〜第６Ｄ図にプロットしてみる
と、少くともグラフ上はｐ１〜ｐ８に真値を得た結果と
、完全に一致してみえる。従ってこのｆｉｔｔｉｎｇは
十分にすぐれたものである。

次に、トランジスタ等価回路定数決定問題のための本装
置の構成例を第８図に示す。以下これについて概要を説
明する。大きな正方形１〜８で表わしたのが、８個の回
路定数ｐ、〜ｐ８に対応する論理セルであり、一般に非
線形な入出力応答Ｖ−ｇ　（ｖ）をもつものとする。各
論理セルの入力端子には、それぞれアンプＡＩ、Ａ２．
・・・Ａ８と抵抗Ｒ１，Ｒ２，・・・Ｒ３、容量Ｃ１，
Ｃ２，・・・Ｃ８が接続され、アンプＡＩ、Ａ２．・・
・Ａ８はλ１．λ２．・・・Ａ８に相当する可変利得を
もつ。

既に述べた通り、このλを調節し適正化することが、本
装置の計算を成功させる重要ポイントである。

容ｆｆ１ｃ１．Ｃ２，・・・Ｃ８は、時間微分項に対応
する電流Ｃ（ｄｖ／ｄｔ）を与える。抵抗は緩和時定数
τ−ＣＲを与える。アンプＡｌ、Ａ２．・・・Ａ８の入
力端子に接続する電流源Ｓｌ、Ｓ２．・・・Ｓ８は、後
述のごとく、式（３１−３）の定数項す、に対応する。

各論理セルの出力に接続する８個のボックス１１゜２１
、・・・８１；・・・・・・１８．２８．・・・８８は
Ｔ　ＩＩＶ　ｒなる電流源を表わし、第ｊパラメータｐ
、を規格化単位ｒ。

で規格化した量ｑ、を、第ｊセルの出力電圧Ｖ。

と対応させる。以上の回路構成により、第８図の回路網
の過渡応答をきめる式はつぎのように書ける。

Ｃ＋　　（ｄｖ＋　　／ｄ　　ｔ）ｗ−Ｖｌ　　／Ｒ＋
ＶＩ　−ｇ　　（Ｖ＋　　）　　　　　　　　　　　−
（４０）ただしｉは論理セルの番号、ｖ、、−Ｖｌは各
々ｉ番目のセルの人、出力電圧を表わす。ＲＣ２ｌｌは
各々抵抗、容量、電流源を表わす。

Ｔ１は第ｉセルの入力と第ｊセルの出力の間の結合強度
を表わし、関数ｇは各セルの人出力応答関数である。式
（３９）　、　（４０）において、つぎの条件を与える
と、 ｄｑ／ｄｔ−−λ（Ａｑ−ｂ） 式に帰着することが分る。

（１）セル間結合強度Ｔ、を容量　Ｃ＋で割った値をａ
ｌｌと対応させる。

（２）電流Ｉ＋を容量Ｃ１で割った値を定数す。

と対応させる。

（３）入出力応答関数ｇを、単に線形関数Ｖ、−ｖ、で
置き換える。

（４）時定数Ｃ＋Ｒ＋は、式 ％式％ ので、Ｒ５−ψとする。

すなわち、以上（１）〜（４）の対応づけにより、第８
図の回路網を用いて式 ｄ　ｑ／ｄ　ｔ−一λ（Ａｑ−ｂ）の計算が自動的に行
われる訳である。

但し、上記以外に、実際には電圧依存性をもつ電流ｇＴ
ｚＶｊおよびＩ、をアナログ的に実現するため、第８図
の右端に示すような制御装置が必要である。この装置は
、まずｒｌｔｔｌｎｇのための参照データとして、ｙパ
ラメータの測定装置により提供される測定データ１０４
を必要とする。次にこの測定データと、各論理セルの出
力電圧とを用いて、式（３１−１）、（３１−２）、（
３１−３）に従い、微分演算、かけ算、たし算を行う演
算装置１０３が必要である。これらの計算はデジタル、
アナログいずれの方式で行なってもよいが、計算精度を
考慮すればデジタル方式が好ましいと考えられる。

以上により作成された式（３１−１）、（３１−２）、
（３１−３）の諸データがデジタル量である場合は、こ
れをＤ−Ａ変換装置１０１を通してアナログ信号に変換
し、さらにアドレスユニット１００を用いて所定の番地
（１，ｊ）・・・・・・・・・ｉ、ｊとも１から８まで
変化する・・・・・・・・・にこれを伝達する。

以上がトランジスタ等価回路定数決定問題を解くための
、本発明の主旨にもとづく装置の一例である。

［発明の効果］ 以上述べたように本発明の方法によれば、行列方程式Ａ
ｘｍｂの解を求める際、感度係数λを導入して、これを
微分方程式 ｄ　ｘ　／　ｄ　ｔ−−λ（Ａｘ−ｂ）の積分問題に置
き換えることにより、逆行列計算ｘ−Ａ−’ｂを行うこ
となく、その目的を達することができる。このことは、
行列Ａが大形サイズで、かつゼロ要素が多い場合、とく
に有効な解法となる。

実際の物理系の挙動を記述する微分方程式から派生する
行列方程式の大多数は、このような大形かつゼロ要素が
多いという特徴を備えているから、本解法がこれらに対
して有効に適用されることは、いうまでもない。

本発明の装置によれば、行列方程式Ａｘ−ｂの解の計算
が、微分方程式 ｄ　ｘ　／　ｄ　ｔ　＝−λ（Ａｘ−ｂ）の積分問題に
おきかえられ、電気回路で模倣した装置により実行され
る。これによって回路が定常状態に至れば自然かつ速や
かに解が求められる。本装置は、特に、行列Ａのサイズ
が大きくかつゼロ要素が多い場合（実際の物理系の挙動
を既述する微分方程式から派生する行列は大多数この種
のものである）に特に有効な解法を与える。

【図面の簡単な説明】

第１図は本発明の第１実施例に係る半導体デバイスの動
作解析法を示すフローチャート、第２図は本発明の第２
実施例に係る半導体デバイスの動作解析法を示すフロー
チャート、第３図は本発明の第１実施例に係る半導体デ
バイスの動作解析装置を示す回路図、第４Ａ図は半導体
デバイスの等価回路定数決定問題のための回路図、第４
Ｂ図は第４Ａ図の等価回路図、第５図は本発明の第３実
施例に係る半導体デバイスの等価回路定数決定問題の解
析法を示すフローチャート、第６Ａ図ないし第６Ｄ図は
真値と試行値に対応する４個のｙパラメータを複素平面
上にプロットしたグラフ、第７図は本発明の第３実施例
に係る半導体デバイスの等価回路定数決定問題の解析法
を実行したデータを示す図、第８図は本発明の第２実施
例に係る半導体デバイスの等価回路定数決定問題の解析
装置を示す回路図、第９図は半導体デバイス空間のメツ
シュ分割を示す図、第１０図は電子１、正孔の暢送方程
式及びポアッソン方程式からなる連立方程式系の行列・
ベクトル方程式を示す図である。 １．２．・・・１２・・・論理セル、１１．＋２．・・
・＋２１２・・・電流検出手段、１０（１・・・アドレ
スユニット、１［１１・・・Ｄ／Ａコンバータ、１０３
，１０４・・・回路。

Claims (8)

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