JPH08335225A

JPH08335225A - 半導体デバイスシミュレーション方法および装置

Info

Publication number: JPH08335225A
Application number: JP7164679A
Authority: JP
Inventors: Koichi Fukuda; 浩一福田
Original assignee: Oki Electric Industry Co Ltd
Current assignee: Oki Electric Industry Co Ltd
Priority date: 1995-06-07
Filing date: 1995-06-07
Publication date: 1996-12-17
Anticipated expiration: 2019-09-29
Also published as: JP3569750B2

Abstract

(57)【要約】【構成】電界を求めるポアソン方程式を解き、その解
に基づきボルツマン方程式を解いてキャリアの分布関数
Ｆを求め、この分布関数の値（Ｆ）から求められるキャ
リア濃度をポアソン方程式に代入し、これら両方定式の
演算を繰り返により、集束する分布関数Ｆによって半導
体デバイスの電気特性を求める方法であって、ボルツマ
ン方程式の解法にセルオートマトン法を利用する。【効果】ボルツマン方程式の解法に、初期値等の諸条
件は仮定的に導入されるが、反復演算処理により、これ
ら諸条件は、適性値に向けて順次更新されることから、
仮定についての妥当性が問題となることはなく、また、
演算処理に多大の時間を要することはなく、確率の小さ
な現象による解のばらつきが防止されることから、現実
に即した安定した適正なシミュレーション結果を迅速に
得ることができる。

Description

【発明の詳細な説明】

【０００１】

【産業上の利用分野】本発明は、半導体デバイスの設計
等に使用するのに好適な半導体デバイスシミュレーショ
ン方法および装置に関する。

【０００２】

【従来の技術】半導体デバイスの設計等では、模擬的に
半導体デバイスの電気特性を求めるための半導体デバイ
スシミュレーション装置が用いられている。半導体デバ
イスの電気特性を模擬的に求めるには、不純物分布およ
び境界条件等に沿ってポアソン方程式を解いて電位を求
め、この電位から導き出される電界に基づき、ボルツマ
ン方程式を解き、これにより求められるキャリアの分布
関数が集束するまで両演算が反復される。この分布関数
の集束により、キャリア分布関数を求めることができ、
キャリア分布関数の値である分布関数値から、半導体デ
バイスの電気特性の基になるキャリア濃度分布を求める
ことができ、またこのキャリア濃度分布から、半導体デ
バイスの電圧−電流特性を始めとする、種々の電気特性
を得ることができる。

【０００３】ところで、このような方程式による演算処
理にはコンピュータが利用されており、従来のシミュレ
ーション装置にはポアソン方程式を解く演算処理部およ
びボルツマン方程式を解く演算処理部が設けられてい
る。両方程式は偏微分方程式で表わされ、特にボルツマ
ン方程式は、時間、位置ベクトルおよび波数ベクトルの
合計７つの変数の偏微分方程式で表わされることから、
このボルツマン方程式を比較的容易に解くために、流体
モデル法およびモンテカルロ法を利用することが提案さ
れている（「プロセス・デバイスシミュレーション技
術」、ｐ９１〜１３３およびｐ１９７〜２１４、檀良
編著、産業図書株式会社発行、昭和６３年）。

【０００４】

【発明が解決しようとする課題】しかしながら、流体モ
デルによってボルツマン方程式を解く方法では、例えば
キャリアの温度を仮定し、あるいはキャリアのエネルギ
ー分布関数がマクスウェル−ボルツマンの分布関数であ
る等、方程式の簡略化のために多数の仮定を導入する必
要がある。そのため、この多数の仮定の導入により、現
実の環境での結果と不一致を生じることがあり、常に導
入した仮定についての妥当性が問題となる。また、モン
テカルロ法によってボルツマン方程式を解く方法では、
方程式の簡略化およびそれに伴う多数の仮定の導入が行
なわれることはないが、解法自体が乱数に基づくことか
ら、確率の小さな現象の演算にも時間を要し、演算処理
に多大の時間が掛かる。しかも、確率の小さな現象によ
り、解そのものが引きずられ、このために現実には有り
得ないような値を含む結果が求められることがあり、常
に安定した結果を求めることが困難である。

【０００５】そこで、半導体デバイスの電気特性を模擬
的に、比較的速く、しかも現実に即する安定した結果と
して得ることのできる半導体デバイスシミュレーション
方法および装置の現出が望まれていた。

【０００６】

【課題を解決するための手段】本発明は、基本的には、
ボルツマン方程式を解くにあたり、いわゆるセルオート
マトン（cellular automaton, 自己増殖機械）法を利用
する。先ず、本発明の基本概念を述べる。ボルツマン方
程式では、分布関数が時間、位置ベクトルおよび波数ベ
クトルの関数として表わされる。

【０００７】これに対応して、本発明では、位置ベクト
ルに係わる位置空間および波数ベクトルに係わる運動量
空間をそれぞれ生成する。位置空間を小領域に区画し、
各小領域にそれぞれの小領域におけるキャリアの存在確
率を示すキャリアの分布関数値を代表する代表点を設け
る。また、位置空間の各代表点毎にそれぞれ運動量空間
を対応させ、各運動量空間を小領域に区画し、運動量空
間の各小領域にそれぞれの小領域におけるキャリアの存
在確率を示すキャリアの分布関数値を代表する代表点を
設ける。各位置空間の代表点に初期値を与え、この初期
値に基づき対応する運動量空間に初期値を与える。ま
た、この位置空間の初期値に対応してキャリアを各運動
量空間の各小領域に配分する。

【０００８】各運動量空間でキャリアが散乱されるまで
の時間Δｔを求め、後述するボルツマン方程式の衝突ま
での走行に関する項についての解を求めるために、運動
量空間の時間Δｔ後の運動量空間での運動量の変化量を
求め、さらに、この運動量空間に対応する位置空間での
位置の変化量を求める。これら両変化量から、位置空間
および運動量空間における各小領域でのキャリア存在確
率すなわち各代表点における分布関数値が補正され、時
間Δｔ後の各小領域でのキャリア存在確率が求められ
る。

【０００９】また、ボルツマン方程式の散乱に関する項
についての解を求めるために、運動量空間においてキャ
リアの散乱後の運動量を求め、この運動量空間における
散乱後の運動量の変化によるキャリア存在確率である分
布関数値を求め、この運動量空間の分布関数値から対応
する位置空間での分布関数値を求める。さらに、この位
置空間および運動量空間の各小領域でのキャリア存在確
率を示す各分布関数値を初期値として、再び、散乱され
るまでの時間Δｔを求め、時間Δｔ後の走行に関する項
についての解および散乱に関する項についての解を求め
るステップが所定回数繰り返される。以上がボルツマン
方程式を解くにあたっての本発明の方法および装置につ
いての基本概念である。

【００１０】本発明に係る半導体デバイスシミュレーシ
ョン方法は、上記の基本概念に沿って半導体デバイスの
電気的特性をコンピュータを用いて模擬的に求める半導
体デバイスシミュレーション方法において、対象とする
半導体デバイスの不純物濃度分布および境界条件につい
ての条件を入力し、入力条件に沿って、電位とキャリア
濃度との関係を示すポアソン方程式を解いて電界に関す
る情報を求め、ポアソン方程式を解くことにより求めら
れた電界に基づき、ボルツマン方程式を解いてキャリア
の分布関数を求め、ボルツマン方程式を解くことにより
求められたキャリアの分布関数値から求められるキャリ
ア濃度をポアソン方程式に代入して両方定式の演算を繰
り返し、集束する分布関数に基づいて対象とする半導体
デバイスの電気特性を求める。

【００１１】ボルツマン方程式の解は、以下のステップ１．メッシュ状に区画された多数の小領域からなり、キ
ャリア分布関数を求めようとする材料の形状に対応した
位置空間を生成し、該位置空間の各小領域にそれぞれの
領域の分布関数値（Ｆ）を代表する代表点を設け、該代
表点のそれぞれに分布関数値（Ｆ）の初期値を与え、２．位置空間の各代表点毎にメッシュ状に区画された多
数の小領域からなる運動量空間を生成し、該運動量空間
の各小領域に代表点を設け、該代表点のそれぞれに、位
置空間の各代表点に与えられた分布関数値（Ｆ）の初期
値に基づく分布関数値（ｆ）の初期値を与え、３．運動量空間の分布関数値（ｆ）に対応してキャリア
を対応する運動量空間に配分し、４．該運動量空間に配分されたキャリアが散乱されるま
での、半導体材料のエネルギーバンド構造に応じて決ま
る時間Δｔを求め、５．ポアソン方程式の解から求められた電界に基づい
て、エネルギーバンド構造により全運動量空間のキャリ
アの加速運動による時間Δｔ後の位置および運動量につ
いての各変化量を求め、６．キャリアの運動後の位置空間および運動量空間上の
位置を基に、各空間の各小領域の前記代表点における分
布関数値（Ｆ、ｆ）を補正し、７．運動量空間において、散乱の種類に応じたキャリア
散乱後の運動量を求め、当該運動量空間上の位置を基
に、該運動量空間の各小領域の各代表点における分布関
数値（ｆ）を補正し、８．分布関数値（Ｆ、ｆ）が補正された位置空間および
運動量空間について、前記３．から前記７．のステップ
を所定回数繰り返し、分布関数値（Ｆ）から求められ、
ポアソン方程式に代入されるキャリア濃度を逐次更新す
る、の実行により求められることを特徴とする。

【００１２】また、本発明に係る半導体デバイスシミュ
レーション装置は、対象とする半導体デバイスの不純物
濃度分布および境界条件を含む条件が入力される入力装
置と、該入力装置からの入力条件に基づいてキャリア分
布関数についての演算を行うプロセッサと、該プロセッ
サから出力されるキャリア分布関するに情報を表示する
出力装置とを含む。

【００１３】プロセッサは、入力装置からの条件に沿っ
て電界に関するポアソン方程式を解く演算処理部と、ポ
アソン方程式を解くことにより求められた電界に基づ
き、キャリアの分布関数を求めるボルツマン方程式を解
く演算処理部と、両演算処理部の反復演算処理により求
められたキャリア分布関数が集束するか否かを判定する
集束判定部とを備える。

【００１４】ボルツマン方程式を解く演算処理部は、メ
ッシュ状に区画された小領域からなる位置空間を生成す
る位置空間生成部分と、位置空間の各小領域に対応する
複数の運動量空間であってそれぞれがメッシュ状に区画
された小領域からなる運動量空間を生成する運動量空間
生成部分と、位置空間の各小領域に与えられたキャリア
の分布関数値（Ｆ）に基づく分布関数値（ｆ）が与えら
れた運動量空間に位置空間での分布関数値（Ｆ）に対応
して配分されたキャリアが散乱されるまでの時間Δｔを
求めるモデル演算処理部分と、該演算処理部分により求
められた時間Δｔ後のキャリアの位置空間および運動量
空間における位置および運動量についての変化量を求
め、該変化量に基づいて位置空間および運動空間の各小
領域の分布関数値（Ｆ、ｆ）を補正する走行項演算処理
部分と、位置空間の各小領域の補正された分布関数値
（Ｆ）に対応する運動量空間において、散乱の種類に応
じたキャリア散乱後の運動量を求め、当該運動量空間上
の位置を基に、該運動量空間の各小領域の分布関数値
（ｆ）を補正する散乱項演算処理部分とを有する。

【００１５】集束判定部は、走行項演算処理部分および
散乱項演算処理部分により補正された位置空間および運
動量空間における各小領域での分布関数値（Ｆ、ｆ）を
もとに、モデル演算処理部分、走行項演算処理部分およ
び散乱項演算処理部分の所定回数の反復演算処理および
ポアソン方程式を解く演算処理部への反復演算処理を経
て得られる分布関数値（Ｆ）から求められるキャリアの
分布関数が集束するか否かを判定する。

【００１６】

【作用】本発明の方法では、位置空間の各代表点に与え
られる初期値は、従来の流体モデル法におけると同様
に、仮定的に導入されるが、この流体モデル法における
ように仮定条件が固定的に取り扱われることはなく、ボ
ルツマン方程式の反復演算およびポアソン方程式の演算
をも含む反復演算によって、これら仮定的に導入される
条件等は、適性値に向けて順次更新される。従って、従
来のように特定の仮定についての妥当性が問題となるこ
とはない。また、乱数に基づく演算処理ではないことか
ら、演算処理にモンテカルロ法におけるような多大の演
算処理時間を要することはなく、しかも、確率の小さな
現象による解のばらつきが防止される。

【００１７】また、本発明の装置では、該装置を構成す
る各部をコンピュータにより構成することにより、本発
明の方法を適切かつ有効に実施することができる。これ
により、従来に比較して高い精度でしかも高速に適正な
半導体デバイス特性のシミュレーション結果を得ること
ができる。

【００１８】

【実施例】以下、本発明を図示の実施例に沿って詳細に
説明する。図１は、本発明に係るシミュレーション装置
のブロック図である。本発明に係る半導体デバイスシミ
ュレーション装置１０は、対象とする半導体デバイスの
構造的、電気的諸条件および必要とされる初期条件等が
入力される入力装置１１と、該入力装置からの情報に基
づいて半導体デバイスのキャリア分布関数についての演
算を行うプロセッサ１２と、該プロセッサから出力され
るキャリア分布関数に関する情報を表示する出力装置１
３とを含み、プロセッサ１２には、その演算処理のため
に必要な情報を格納するメモリ１４が接続されている。

【００１９】入力装置１１から入力される構造的あるい
は電気的条件としては、対象とする半導体デバイスが例
えばＭＯＳトランジスタのような素子であるとき、ゲー
ト長、ゲート幅、ゲート絶縁膜の厚さのような形状に関
する条件、半導体基板にソースあるいはドレイン領域を
形成するためにこの基板に導入される不純物の濃度分
布、接合部における境界条件等の解析に必要な条件があ
る。また初期条件として、キャリア濃度（ｐ−ｎ：正孔
と電子の差）が不純物濃度（Ｎ_D −Ｎ_A ：ドナー濃度と
アクセプタ濃度との差）に等しいと仮定する電荷中性条
件のような初期値条件が入力される。

【００２０】プロセッサ１２は、入力装置１１からの条
件に沿って、電位とキャリア濃度の関係を示すポアソン
方程式（１式）を解く演算処理部１５と、ポアソン方程
式を解くことにより求められた電界に基づき、キャリア
の分布関数を求めるボルツマン方程式（２式）を解く演
算処理部１６と、両演算処理部１５、１６の反復演算処
理により求められたキャリア分布関数が集束するか否か
を判定する集束判定部１７とを備える。

【００２１】ポアソン方程式の演算処理部１５は、例え
ば初期条件として入力された前記電荷中性条件に沿った
電位を算出する。この電位は、引き続くボルツマン方程
式の演算処理部１６のために微分により、電界に変換さ
れる。

【００２２】この電界に基づいてキャリア分布関数を求
める演算処理部１６は、メッシュ状に区画された小領域
からなる位置空間を生成する位置空間生成部分１８と、
位置空間の各小領域に対応する複数の運動量空間であっ
てそれぞれがメッシュ状に区画された小領域からなる運
動量空間を生成する運動量空間生成部分１９と、運動量
空間でキャリアが散乱されるまでの時間Δｔを求めるモ
デル演算処理部分２０と、ボルツマン方程式（２式）の
左辺で示されるキャリア走行項についての演算処理を行
う走行項演算処理部分２１と、ボルツマン方程式（２
式）の右辺で示されるキャリア散乱項についてキャリア
の散乱後の運動量を求める散乱項演算処理部分２２とを
有する。

【００２３】ボルツマン方程式（２式）に示されるキャ
リア分布関数Ｆは、位置についての位置ベクトルおよび
運動量についての波数ベクトルの関数であり、これに関
連して、位置空間生成部分１８および運動量空間生成部
分１９は、それぞれ位置空間および運動量空間をメモリ
１４に生成する。

【００２４】図２は位置空間を概念的に示す説明図であ
り、図３および図４はそれぞれ運動量空間を概念的に示
す説明図である。図示の例では、位置空間２３および運
動量空間２４は、説明の簡略化のために、それぞれ直交
するメッシュによって多数の小領域に区画された２次元
メッシュで表現されている。これら空間２３および２４
を３次元で表わすことができる。また、位置空間２３の
実線は、キャリア分布関数を求めようとする材料の形状
に対応している。

【００２５】位置空間２３の各小領域は、キャリアの存
在確率である分布関数値（Ｆ）をそれぞれ与えられ、各
小領域の分布関数値（Ｆ）を代表する代表点が例えばメ
ッシュの交点あるいは小領域の中央部等にそれぞれ設定
される。図２には、図面の簡略化のために小領域２３Ａ
の分布関数値（Ｆ）を代表する１つのキャリアＣ１のみ
が示されている。位置空間２３の１つのキャリアＣ１に
対して、例えば図３で示したような１つの運動量空間２
４が対応する。従って、位置空間２３の小領域の数に対
応した数の運動量空間２４が生成される。

【００２６】各運動空間２４の小領域には、該運動空間
に対応する位置空間２３のキャリアＣ１の分布関数値
（Ｆ）に対応した分布関数値（ｆ）が与えられる。例え
ば、キャリアＣ１の分布関数値（Ｆ）が１であったとす
ると、このキャリアＣ１に対応する運動量空間２４の全
小領域における分布関数値（ｆ）の和すなわち対応する
運動量空間２４の全キャリアＣ２の分布関数値（ｆ）の
和が、対応するキャリアＣ１の分布関数値（Ｆ）である
１に等しくなるように設定される。

【００２７】この運動量空間２４のキャリアＣ２の数に
ついては、複数のキャリアＣ２にそれぞれ等価の重みを
与えるか、あるいは単一のキャリアＣ２にこれに対応す
るキャリアＣ１の分布関数値（Ｆ）に応じた分布関数値
（ｆ）を与えることができる。複数のキャリアＣ２に等
価の重みを与える場合、例えばキャリアＣ１の分布関数
値（Ｆ）が１であり、これに対応するキャリアＣ２が図
３に示されているように２個であるとき、各キャリアＣ
２には、それぞれ０．５の重みが与えられる。従って、
０．５の重みが与えられたキャリアＣ２は、これに与え
られた分布関数値（ｆ）の半値の存在確率を持つことに
なる。

【００２８】図４には運動量空間２４の単一のキャリア
Ｃ２が示されており、この単一のキャリアＣ２には１．
０の重みが与えられている。従って、この単一のキャリ
アＣ２は、対応するキャリアＣ１の分布関数値（Ｆ）に
応じて与えられた分布関数値（ｆ）を持つこととなる。
次に述べるキャリア走行項についての運動量空間２４で
各キャリアの運動を追う必要があることから、このよう
に各運動量空間２４毎に単一のキャリアＣ２を適用する
ことが、精度の低下を招くことなく演算処理時間の短縮
化を図る上で、有利である。

【００２９】ボルツマン方程式（２式）の左辺で示され
るキャリア走行項は、キャリアが衝突すなわち散乱する
までの時間Δｔの間におけるキャリアの走行運動を表わ
す項である。この時間Δｔを求めるモデル演算処理部分
２０は、半導体材料に応じたキャリアの物理モデルに応
じて、例えば格子散乱、不純物散乱等毎に散乱が生じる
までの時間Δｔを算出する。

【００３０】また、ボルツマン方程式（２式）の右辺で
ある散乱項を零におき、この式に、ポアソン方程式演算
処理部１５で求められた電界およびモデル演算処理部分
２０で求められた時間Δｔを代入して演算処理すること
により、散乱を生じるまでのキャリアの運動量空間２４
での運動を求めることができる。走行項演算処理部分２
１は、キャリアに応じたエネルギバンドモデルを用い
て、運動量空間２４でのキャリアＣ２の運動量の変化量
を求め、また、求めた運動量空間２４での変化量をもと
に位置空間２３での対応するキャリアの位置についての
変化量を求める。このキャリアの位置の変化量は、エネ
ルギーバンドモデルと電界とから求めることもできる。

【００３１】図４には、運動量空間２４で黒丸位置から
白丸位置に移動したキャリアＣ２が示されている。走行
項演算処理部分２１は、このキャリアＣ２の変化量を算
出すると共に、この移動によって生じる運動量空間２４
の各小領域での分布関数値（ｆ）を補正する。

【００３２】この補正については、図４を参照するに、
キャリアＣ２の分布関数値（ｆ）は、移動後のキャリア
Ｃ２を取り巻く４つの代表点Ａ、Ｂ、Ｃ、Ｄおよびキャ
リアＣ２の移動後の点で区画される面積Ｓ１、Ｓ２、Ｓ
３、Ｓ４に応じて、次式Ａ＝Ｓ３／ＳＢ＝Ｓ４／ＳＣ＝Ｓ１／ＳＤ＝Ｓ２／Ｓ（各式において、ＳはＳ１，Ｓ２，Ｓ３，Ｓ４の和であ
る。）で示される割合で、近接する各代表点Ａ、Ｂ、
Ｃ、Ｄに分け与えられる。従って、走行項演算処理部分
２１は、各代表点をそれぞれの新たな分布関数値（ｆ）
に補正する。なお、分布関数値（ｆ）は、図４に沿って
説明したこの補正に先立って、これと同様の方法によ
り、キャリアの位置空間での移動に応じて、走行項演算
処理部分２１により予め補正されている。

【００３３】位置空間２３あるいは運動量空間２４が１
次元で表わされている場合、分布関数値（Ｆ、ｆ）の補
正は、先に述べた面積比におけると同様に、近接する代
表点間での距離比に応じて補正され、また位置空間２３
あるいは運動量空間２４が３次元で表わされている場
合、分布関数値（Ｆ、ｆ）の補正は、体積比に応じて同
様に補正される。また、図３に示したように、複数の等
価のキャリアＣ２が適用されている場合、走行項演算処
理部分２１は、各キャリアＣ２が移動先で取り囲まれ
た、それぞれ斜線で示す小領域の代表点に、分布関数値
（ｆ）をキャリアの重みで等分した値を与える。

【００３４】ボルツマン方程式（２式）の右辺で示され
るキャリア散乱項についてキャリアの散乱後の運動量を
求める散乱項演算処理部分２２は、運動量空間２４で生
じる各種の散乱についてその散乱後のキャリアの運動量
を求め、この運動量空間２４での運動量に基づいて、走
行項演算処理部分２１におけると同様に、運動量空間２
４の各小領域での代表点の分布関数値（ｆ）を補正す
る。キャリアの散乱後の運動量を求めるについて、散乱
の種類を確率分布に基づく乱数に応じて選択することが
できる。しかしながら、より迅速かつ高精度の演算処理
を求める上で、各散乱モデルについての確率を求め、こ
の確率に基づいてキャリア運動量を算出することが望ま
しい。

【００３５】図５は、運動量空間２４におけるキャリア
Ｃ２の散乱での運動量変化による分布関数値（ｆ）の補
正についての説明図である。黒丸で示されるキャリアＣ
２が確率Ｐ１で生じる散乱により白丸の位置（Ｃ２−
１）に移動し、また確率Ｐ２で生じる他の散乱により２
重の白丸の位置（Ｃ２−２）に移動したとする。このと
き、キャリアＣ２−１には、キャリアＣ２が有していた
分布関数値（ｆ）に確率Ｐ１の重みを掛けた値が与えら
れ、またキャリアＣ２−２には、キャリアＣ２が有して
いた分布関数値（ｆ）に確率Ｐ２の重みを掛けた値が与
えられる。

【００３６】このキャリアの散乱による補正は、図４に
ついて説明したと同様に、キャリアＣ２−１について
は、キャリアＣ２の分布関数値（ｆ）に確率Ｐ１を掛け
た値が、キャリアＣ２−１を取り巻く４つの代表点Ａ、
Ｂ、Ｅ、Ｆにそれぞれ区画面積Ｓ１、Ｓ２、Ｓ３、Ｓ４
に応じて、次式Ａ＝Ｓ４／ＳＢ＝Ｓ３／ＳＥ＝Ｓ２／ＳＦ＝Ｓ１／Ｓ（各式において、ＳはＳ１，Ｓ２，Ｓ３，Ｓ４の和であ
る。）で示される割合で、近接する各代表点Ａ、Ｂ、
Ｅ、Ｆに分け与えられる。

【００３７】また、キャリアＣ２−２についても、キャ
リアＣ２の分布関数値（ｆ）に確率Ｐ２を掛けた値が、
キャリアＣ２−２を取り巻く４つの代表点Ｂ、Ｃ、Ｄ、
Ｅに同様な区画面積比（Ｓ１，Ｓ２，Ｓ３，Ｓ４）に応
じて分け与えられる。これらの分布関数値（ｆ）の配分
により、運動量空間２４の分布関数値（ｆ）の散乱によ
る補正がなされ、この運動量空間２４における各代表点
の分布関数値の補正により、ボルツマン方程式の基本的
な１サイクルの演算処理が終る。

【００３８】この１サイクルの演算処理後、走行項演算
処理部分２１により求められた位置空間２３における各
分布関数値および散乱項演算処理部分２２により求めら
れた運動量空間２４における各分布関数値を、それぞれ
両空間２３および２４の初期値として、モデル演算処理
部分２０、走行項演算処理部分２１および散乱項演算処
理部分２２での前記したと同様な演算処理サイクルが所
定回数繰り返される。

【００３９】モデル演算処理部分２０、走行項演算処理
部分２１および散乱項演算処理部分２２の反復演算が所
定回数に達すると、集束判定部１７が走行項演算処理部
分２１により求められた位置空間２３における各分布関
数値（Ｆ）から求められるキャリアの分布関数Ｆが集束
したか否か、すなわち分布関数Ｆがほぼ誤差内で集束し
たか否かを判定する。集束していなければ、その分布関
数値（Ｆ）で求められるキャリア濃度を演算処理部１５
におけるキャリア濃度として、再度、各演算処理部１
５、１６での演算処理が繰り返される。

【００４０】出力装置１３は、この集束するキャリア分
布関数Ｆを解として表示する。この出力装置１３に、キ
ャリア分布関数Ｆから求められる種々の電気特性に関す
る情報を必要に応じて表示させることができる。

【００４１】図６および図７は、本発明のシミュレーシ
ョン方法の手順を示すフローチャートである。半導体デ
バイスシミュレーション装置１０の作動の概略を図６に
沿って説明する。

【００４２】先ず、入力装置１１に、対象とする半導体
デバイスの不純物濃度分布、境界条件あるいは電荷中性
条件等の初期条件が設定値設定として入力される（ステ
ップＳ６１）。これらの入力条件に沿って、ポアソン方
程式を解く演算処理部１５が演算処理を行い、電界に関
する情報を求める（ステップＳ６２）。続いて、ポアソ
ン方程式を解く演算処理部１５により求められた電界に
基づき、ボルツマン方程式を解く演算処理部１６がセル
オートマトン法に沿って、キャリア分布関数Ｆを求める
（ステップＳ６３）。演算処理部１６により求められた
キャリア分布関数Ｆが、集束するか否かを集束判定部１
７により、判定される（ステップＳ６４）。

【００４３】キャリア分布関数Ｆが集束していない場
合、このキャリア分布関数Ｆから求められるキャリア濃
度が演算処理部１５のキャリア濃度として、再度ステッ
プＳ６２以下の演算処理が繰り返される。また、キャリ
ア分布関数Ｆが集束する場合、出力装置１３により、キ
ャリア分布関数Ｆに関する情報が出力される。

【００４４】先に述べたボルツマン方程式を解く演算処
理（ステップＳ６３）では、図７に示されているよう
に、位置空間生成部分１８および運動量空間生成部分１
９によって位置空間および運動空間がそれぞれ生成され
かつそれぞれがメッシュ状に区画され、各空間の小領域
について代表点が決められ、位置空間の個々の代表点に
ついての分布関数値（Ｆ）が初期値として設定され、こ
れに基づいて運動量空間の各初期値が決まる（ステップ
Ｓ７１）。

【００４５】この位置空間の初期値に、マクスウェル−
ボルツマン分布関数を適用することができ、また、従来
の前記した流体モデル法によって得られた分布関数を適
用することができる。この流体モデル法によって得られ
た分布関数は、適正な分布関数に近似する関数であるこ
とから、これに基づき各代表点の初期値を設定すること
により、すなわち、セルオートマトン法の初期解に流体
モデル法を用いることにより、演算処理時間の短縮化を
図ることができる。

【００４６】ステップＳ７１に示される初期設定は、図
７に示されたサブルーチンを実行する際の最初に実行さ
れるのみであり、集束ステップＳ６４からの繰り返しに
よっる二回目以降のループではこのステップＳ７１が繰
り返されることはない。

【００４７】また、運動量空間生成部分１９は、運動量
空間の分布関数値（ｆ）に対応して、キャリアを対応す
る運動量空間２４に配分する（ステップＳ７２）。運動
量空間２４に配分されたキャリアが散乱するまでの時間
Δｔがモデル演算処理部分２０により求められ（ステッ
プＳ７３）、この時間Δｔ後のキャリアの位置空間２３
および運動量空間２４での変化量が走行項演算処理部分
２１により求められ（ステップＳ７４）、その結果に基
づいて位置空間２３および運動量空間２４での分布関数
値（Ｆ、ｆ）が補正される（ステップＳ７５）。

【００４８】続いて、散乱項演算処理部分２２が運動量
空間２４でのキャリア散乱後の運動量を求め、この運動
量に基づいて、運動量空間２４の分布関数値（ｆ）を補
正する（ステップＳ７６）。この補正された位置空間２
３および運動量空間２４での分布関数値（Ｆ、ｆ）に基
づいて、ステップＳ７２〜ステップＳ７６が所定回数反
復され（ステップＳ７７）、その後、図６に示した集束
判断ステップＳ６４に移る。

【００４９】本発明の半導体デバイスシミュレーション
装置１０によれば、位置空間の各代表点に与えられる初
期値を決める仮定条件は、固定的に取り扱われることは
なく、ボルツマン方程式の反復演算およびポアソン方程
式の演算をも含む反復演算によって、適性値に向けて順
次更新される。従って、従来のように特定の仮定につい
ての妥当性が問題となることはない。また、乱数に基づ
く演算処理ではないことから、演算処理にモンテカルロ
法におけるような多大の演算処理時間を要することはな
く、しかも、確率の小さな現象による解のばらつきが防
止される。従って、迅速に、しかも高い精度で適正なシ
ミュレーション結果を得ることができる。

【００５０】図８は、半導体デバイスシミュレーション
装置１０によって得られたキャリア濃度分布特性曲線を
示すグラフである。この例では、基板２５としてＰ型シ
リコンを用いたＭＯＳＦＥＴ２６について、ゲート電極
２７の両側にソース／ドレイン領域２８を形成するため
に、例えばＡｓ＋のようなイオンが注入され、ある外部
電圧が電極に印加されたときのキャリア濃度の分布特性
曲線がＭＯＳＦＥＴ２６の縦断面図上に求められた。そ
れぞれ１０¹⁸〜１０²⁰／cm³ の等濃度を示す曲線が滑ら
かな特性曲線で表わされているように、モンテカルロ法
におけるような統計的揺らぎによる特異な点が出現する
ことなく、より正確なキャリア濃度分布特性曲線が迅速
に得られた。

【００５１】図９は、運動量空間２４の他の例を示す説
明図である。運動量空間２４を、図３〜図５に示したよ
うに方形に区画することに代えて、図９に示すように、
円形メッシュを適用することができる。図では説明の簡
略化のために２次元で示されているが、デカルト座標で
は球形メッシュが適用される。

【００５２】図９で、キャリアが散乱前の運動量を示す
位置Ｃ３からある散乱によってＣ３−１に示した位置へ
移動したとすると、この移動による分布関数値（ｆ）の
補正は次のようになされる。移動後の位置Ｃ３−１を取
り囲む代表点Ａ、Ｂ、Ｃ、Ｄに分布関数値（ｆ）が配分
されるが、先ずＸ軸およびＹ軸上への配分比が運動ベク
トルのＸ軸成分およびＹ軸成分比であるｃｏｓθ、ｓｉ
ｎθに沿って配分される。

【００５３】次にＸ軸上に配分された分布関数値（ｆ）
×ｃｏｓθの値が、Ｘ軸上の代表点Ｃ、Ｄ上に、図４に
沿って説明したと同様な距離比に応じて、Ｃ＝ｎ／（ｎ＋ｍ）Ｄ＝ｍ／（ｎ＋ｍ）の割合で配分される。同様に、Ｙ軸上に配分された分布
関数値（ｆ）×ｓｉｎθの値が、Ｙ軸上の代表点Ａ、Ｂ
上に、その距離比に応じて代表点Ｃ、Ｄにおけると同様
に配分される。このＸ軸、Ｙ軸上の代表点の分布関数値
（ｆ）の配分により、運動量空間２４の分布関数値
（ｆ）の散乱による補正がなされる。

【００５４】このように、円形メッシュあるいは球状メ
ッシュを使用することにより、代表点を軸（Ｘ、Ｙ、ま
たはＺ）上に集約することができ、演算の簡素化および
それによる高速化を図ることができる。

【００５５】また、以上に説明したところでは、各位置
空間および運動量空間のメッシュで区画された各小領域
はその領域内では均等な分布関数値すなわち存在確率を
示すと考えたが、このような仮定で、演算時間の短縮化
を図るためにメッシュを粗く設定すると、疑似拡散と呼
ばれる現象が生じ、誤差が大きくなる。これを防止する
ために、各小領域の存在確率を示す分布関数値を補間関
数により修正することにより、メッシュを粗く設定して
も疑似拡散を防止し、精度の高い演算をより短時間で行
うことが可能になる。

【００５６】図１０は、この補間関数の原理の説明図で
ある。横軸には、各小領域を区画するメッシュのＸ軸が
示され、縦軸には各小領域の代表点が表わす存在確率が
示されている。各小領域の代表点がそれぞれ黒点で示さ
れているが、これら各小領域が均等な存在確率を示すと
仮定すると、分布関数は、図中破線で示すような折れ線
で示されることとなる。しなしながら、このような仮定
は非現実的であり、真の分布関数の近似に過ぎない。そ
こで、各代表点を通る関数曲線を補間関数Ｈとして採用
し、キャリアの分布関数値（Ｆ、ｆ）を隣接する各代表
点に配分する際、この補間関数Ｈにより配分比を修正す
ることにより、メッシュを比較的細かくすることなく疑
似拡散現象を防止することができることから、より短時
間での高精度の演算処置が可能になる。

【００５７】本発明は以上の実施例に限定されない。位
置空間および運動量空間の各Ｘ軸およびＹ軸をそれぞれ
等間隔の等値で区画した例を示したが、各軸の間隔を不
等とすることができる。この場合、円形メッシュは長円
メッシュになる。また、磁場を考慮することができ、こ
の磁場についての演算処理過程をポアソン方程式の演算
処理と、ボルツマン方程式の演算処理との間に、挿入す
ることができる。また、磁場を考慮した場合には、走行
項演算処理部分２１での演算に磁場の効果を考慮する必
要がある。

【００５８】

【発明の効果】以上説明した本発明の方法では、ボルツ
マン方程式を解くに際し、位置空間の各代表点に与えら
れる初期値等は仮定的に導入されるが、ボルツマン方程
式の反復演算およびポアソン方程式の演算をも含む反復
演算によって、これら仮定的に導入される諸条件は、固
定的に取り扱われることはなく適性値に向けて順次更新
される。そのため、従来のような特定の仮定についての
妥当性が問題となることはない。また、多大の演算処理
時間を要することはなく、確率の小さな現象による解の
ばらつきが防止される。

【００５９】従って、本発明の方法によれば、現実に即
する安定した適正なシミュレーション結果が迅速に得ら
れる。また、本発明の装置よれば、該装置を構成する各
部をコンピュータにより構成することにより、本発明の
方法を適切かつ有効に実施することができることから、
適正なシミュレーション結果を、高い精度で、しかも迅
速に得ることができる。

【図面の簡単な説明】

【図１】本発明に係るシミュレーション装置のブロック
図である。

【図２】本発明に係る位置空間の説明図である。

【図３】本発明に係る運動量空間の説明図（その１）で
ある。

【図４】本発明に係る運動量空間の説明図（その２）で
ある。

【図５】本発明に係る運動量空間の説明図（その３）で
ある。

【図６】本発明に係るシミュレーション方法の手順を示
すフローチャート（その１）である。

【図７】本発明に係るシミュレーション方法の手順を示
すフローチャート（その２）である。

【図８】本発明の装置によるシミュレーション結果を示
すキャリア濃度分布特性曲線図である。

【図９】本発明に係る運動量空間の説明図（その４）で
ある。

【図１０】本発明に係る補間関数の説明図である。

【符号の説明】

１０半導体デバイスシミュレーション装置１１入力装置１２プロセッサ１３出力装置１５ポアソン方程式演算処理部１６ボルツマン方程式演算処理部１７集束判定部１８位置空間生成部分１９運動量空間生成部分２０モデル演算処理部分２１走行項演算処理部分２２散乱項演算処理部分２３位置空間２４運動量空間

Claims

【特許請求の範囲】

【請求項１】半導体デバイスの電気的特性をコンピュ
ータを用いて模擬的に求める半導体デバイスシミュレー
ション方法において、対象とする半導体デバイスの不純物濃度分布および境界
条件についての条件を入力し、入力条件に沿って、電位とキャリア濃度との関係を示す
ポアソン方程式を解いて電界に関する情報を求め、ポアソン方程式を解くことにより求められた電界に基づ
き、ボルツマン方程式を解いてキャリアの分布関数を求
め、ボルツマン方程式を解くことにより求められたキャリア
の分布関数値（Ｆ）から求められるキャリア濃度をポア
ソン方程式に代入して両方程式の演算を繰り返し、集束
する分布関数に基づいて対象とする前記半導体デバイス
の電気特性を求める方法であって、ボルツマン方程式の解を以下のステップ１．メッシュ状に区画された多数の小領域からなり、キ
ャリア分布関数を求めようとする材料の形状に対応した
位置空間を生成し、該位置空間の各小領域にそれぞれの
領域の分布関数値（Ｆ）を代表する代表点を設け、該代
表点のそれぞれに分布関数値（Ｆ）の初期値を与え、２．前記位置空間の各代表点毎にメッシュ状に区画され
た多数の小領域からなる運動量空間を生成し、該運動量
空間の各小領域に代表点を設け、該代表点のそれぞれ
に、前記位置空間の各代表点に与えられた分布関数値
（Ｆ）の初期値に基づく分布関数値（ｆ）の初期値を与
え、３．前記運動量空間の分布関数値（ｆ）に対応してキャ
リアを対応する運動量空間に配分し、４．該運動量空間に配分されたキャリアが散乱されるま
での、半導体材料のエネルギーバンド構造に応じて決ま
る時間Δｔを求め、５．ポアソン方程式の解から求められた電界に基づい
て、エネルギーバンド構造により全運動量空間のキャリ
アの加速運動による時間Δｔ後の位置および運動量につ
いての各変化量を求め、６．キャリアの運動後の前記位置空間および前記運動量
空間上の位置を基に、各空間の各小領域の前記代表点に
おける分布関数値（Ｆ、ｆ）を補正し、７．前記運動量空間において、散乱の種類に応じたキャ
リア散乱後の運動量を求め、当該運動量空間上の位置を
基に、該運動量空間の各小領域の前記各代表点における
分布関数値（ｆ）を補正し、８．分布関数値（Ｆ、ｆ）が補正された前記位置空間お
よび前記運動量空間について、前記３．から前記７．の
ステップを所定回数繰り返し、分布関数値（Ｆ）から求
められ、ポアソン方程式に代入されるキャリア濃度を逐
次更新する、の実行により求めることを特徴とする、半
導体デバイスシミュレーション方法。
【請求項２】前記運動量空間を区画するメッシュは、
円形メッシュもしくは球形メッシュである請求項１記載
の半導体デバイスシミュレーション方法。
【請求項３】前記運動量空間の分布関数値（ｆ）に対
応して前記運動量空間に配分される前記キャリアは、各
運動量空間に１つが配分され、各キャリアにはそれぞれ
の小領域の分布関数値（Ｆ）に対応する重みが与えられ
ることを特徴とする請求項１記載の半導体デバイスシミ
ュレーション方法。
【請求項４】散乱の種類に応じたキャリア散乱後の運
動量をそれぞれ求め、散乱の種類に応じたそれぞれの運
動量に基づいて、該運動量空間の各小領域の前記各代表
点における分布関数値（ｆ）を補正することを特徴とす
る請求項１記載の半導体デバイスシミュレーション方
法。
【請求項５】半導体デバイスの電気的特性をコンピュ
ータを用いて模擬的に求める半導体デバイスシミュレー
ション装置であって、対象とする半導体デバイスの不純物濃度分布および境界
条件を含む条件が入力される入力装置と、該入力装置からの入力条件に基づいてキャリア分布関数
についての演算を行うプロセッサと、該プロセッサから出力されるキャリア分布関数に関する
情報を表示する出力装置とを含み、前記プロセッサは、前記入力装置からの条件に沿って電界に関するポアソン
方程式を解く演算処理部と、ポアソン方程式を解くことにより求められた電界に基づ
き、キャリアの分布関数を求めるボルツマン方程式を解
く演算処理部と、両演算処理部の反復演算処理により求められたキャリア
分布関数が集束するか否かを判定する集束判定部とを備
え、前記ボルツマン方程式を解く演算処理部は、メッシュ状に区画された小領域からなる位置空間を生成
する位置空間生成部分と、前記位置空間の各小領域に対応する複数の運動量空間で
あってそれぞれがメッシュ状に区画された小領域からな
る運動量空間を生成する運動量空間生成部分と、前記位置空間の各小領域に与えられたキャリアの分布関
数値（Ｆ）に基づく分布関数値（ｆ）が与えられた前記
運動量空間に前記位置空間での分布関数値（Ｆ）に対応
して配分されたキャリアが散乱されるまでの時間Δｔを
求めるモデル演算処理部分と、該演算処理部分により求められた時間Δｔ後のキャリア
の前記位置空間および前記運動量空間における位置およ
び運動量についての変化量を求め、該変化量に基づいて
前記位置空間および前記運動空間の各小領域の分布関数
値（Ｆ、ｆ）を補正する走行項演算処理部分と、前記運動量空間の各小領域の補正された分布関数値
（ｆ）に対応する前記運動量空間において、散乱の種類
に応じたキャリア散乱後の運動量を求め、当該運動量空
間上の位置を基に、該運動量空間の各小領域の分布関数
値（ｆ）を補正する散乱項演算処理部分とを有し、前記判定部が、前記走行項演算処理部分および前記散乱
項演算処理部分により補正された前記位置空間および前
記運動量空間における各小領域での分布関数値（Ｆ、
ｆ）をもとに、前記モデル演算処理部分、前記走行項演
算処理部分および前記散乱項演算処理部分の所定回数の
反復演算処理および前記ポアソン方程式を解く演算処理
部への反復演算処理を経て得られる分布関数値（Ｆ）か
ら求められるキャリアの分布関数が集束するか否かを判
定することを特徴とする半導体デバイスシミュレーショ
ン装置。