(57)【特許請求の範囲】
【請求項１】
半導体デバイスの電気的特性をコンピュータを用いて模擬的に求める半導体デバイスシミュレーション方法において、
対象とする半導体デバイスの不純物濃度分布、及び電荷中性条件を含む条件を入力し、
入力条件から電位とキャリア濃度との関係を示すポアソン方程式を解いて電界を演算し、
ポアソン方程式を解くことにより求められた電界Ｅに基づき、ボルツマン方程式を解いてキャリアの分布関数を求め、
ボルツマン方程式を解くことにより求められたキャリアの分布関数値（Ｆ）から求められるキャリア濃度を前記ポアソン方程式に代入して両方程式の演算を繰り返し、集束する分布関数に基づいて対象とする前記半導体デバイスの電気特性を求める方法であって、
前記ボルツマン方程式の解として以下のステップ
１．材料の形状に対応してメッシュ状に区画された多数の小領域からなる各位置空間を生成し、該各位置空間に代表点を設け、該各代表点のそれぞれにキャリア存在数を示す初期値を与え、該初期値を前記分布関数Ｆに代入して分布関数値（Ｆ）を算出し、
２．前記位置空間の各代表点毎にメッシュ状に区画された多数の小領域からなる運動量空間を生成し、該各運動量空間に代表点を設け、該代表点のそれぞれにキャリア存在数を示す初期値を与え、該初期値を前記分布関数値（Ｆ）に代入して分布関数値（ｆ）を算出し、
３．前記運動量空間の分布関数値（ｆ）に対応してキャリアを対応する運動量空間に配分し、
４．半導体材料の物理モデルに応じて該運動量空間に配分されたキャリアが散乱されるまでの時間Δｔを求め、
５．ポアソン方程式の解から求められた電界及びエネルギーバンド構造に基づいて全運動量空間のキャリアの加速運動による時間Δｔ後の位置および運動量についての各変化量を求め、
６．キャリアの運動後の前記位置空間および前記運動量空間上の位置を基に、各空間の各小領域の前記代表点における分布関数値（Ｆ、ｆ）を補正し、
７．前記運動量空間において、散乱の種類に応じたキャリア散乱後の運動量を求め、当該運動量空間上の位置を基に、該運動量空間の各小領域の前記各代表点における分布関数値（ｆ）を補正し、
８．分布関数値（Ｆ、ｆ）が補正された前記位置空間および前記運動量空間について、前記３．から前記７．のステップを所定回数繰り返し、分布関数値（Ｆ）から求められてポアソン方程式に代入されるキャリア濃度を逐次更新する、
の実行により求めることを特徴とする、半導体デバイスシミュレーション方法。
【請求項２】
前記運動量空間を区画するメッシュは、円形メッシュもしくは球形メッシュである請求項１記載の半導体デバイスシミュレーション方法。
【請求項３】
前記運動量空間の分布関数値（ｆ）に対応して前記運動量空間に配分される前記キャリアは、各運動量空間に１つが配分され、各キャリアにはそれぞれの小領域の分布関数値（Ｆ）に対応する重みが与えられることを特徴とする請求項１記載の半導体デバイスシミュレーション方法。
【請求項４】
散乱の種類に応じたキャリア散乱後の運動量をそれぞれ求め、散乱の種類に応じたそれぞれの運動量に基づいて、該運動量空間の各小領域の前記各代表点における分布関数値（ｆ）を補正することを特徴とする請求項１記載の半導体デバイスシミュレーション方法。
【請求項５】
半導体デバイスの電気的特性をコンピュータを用いて模擬的に求める半導体デバイスシミュレーション装置であって、
対象とする半導体デバイスの不純物濃度分布および境界条件を含む条件が入力される入力装置と、
該入力装置からの入力条件に基づいてキャリア分布関数についての演算を行うプロセッサと、
該プロセッサから出力されるキャリア分布関数に関する情報を表示する出力装置とを含み、
前記プロセッサは、
前記入力装置からの条件に沿って電位とキャリア濃度に関するポアソン方程式を解く演算処理部と、
ポアソン方程式を解くことにより求められた電界に基づき、キャリアの分布関数を求めるボルツマン方程式を解く演算処理部と、
両演算処理部の反復演算処理によりキャリア分布関数が集束するか否かを判定する集束判定部とを備え、
前記ボルツマン方程式を解く演算処理部は、
メッシュ状に区画された小領域からなる位置空間を生成する位置空間生成部分と、
前記位置空間の各小領域に対応する複数の運動量空間であってそれぞれがメッシュ状に区画された小領域からなる運動量空間を生成する運動量空間生成部分と、
前記位置空間の各小領域に与えられたキャリアの分布関数値（Ｆ）に基づく分布関数値（ｆ）が与えられた前記運動量空間において配分されたキャリアが散乱されるまでの時間Δｔを求めるモデル演算処理部分と、
該演算処理部分により求められた時間Δｔ後のキャリアの前記位置空間および前記運動量空間における位置および運動量についての変化量を求め、該変化量に基づいて前記位置空間および前記運動空間の各小領域の分布関数値（Ｆ、ｆ）を補正する走行項演算処理部分と、
前記運動量空間の各小領域の補正された分布関数値（ｆ）に対応する前記運動量空間において、散乱の種類に応じたキャリア散乱後の運動量を求め、当該運動量空間上の位置を基に、該運動量空間の各小領域の分布関数値（ｆ）を補正する散乱項演算処理部分とを有し、
前記判定部が、前記走行項演算処理部分および前記散乱項演算処理部分により補正された前記位置空間および前記運動量空間における各小領域での分布関数値（Ｆ、ｆ）をもとに、前記モデル演算処理部分、前記走行項演算処理部分および前記散乱項演算処理部分の所定回数の反復演算処理および前記ポアソン方程式を解く演算処理部への反復演算処理を経て得られる分布関数値（Ｆ）から求められるキャリアの分布関数が集束するか否かを判定することを特徴とする半導体デバイスシミュレーション装置。
【発明の詳細な説明】
【０００１】
【産業上の利用分野】
本発明は、半導体デバイスの設計等に使用するのに好適な半導体デバイスシミュレーション方法および装置に関する。
【０００２】
【従来の技術】
半導体デバイスの設計等では、模擬的に半導体デバイスの電気特性を求めるための半導体デバイスシミュレーション装置が用いられている。
半導体デバイスの電気特性を模擬的に求めるには、不純物分布および境界条件等に沿ってポアソン方程式を解いて電位を求め、この電位から導き出される電界に基づき、ボルツマン方程式を解き、これにより求められるキャリアの分布関数が集束するまで両演算が反復される。
この分布関数の集束により、キャリア分布関数を求めることができ、キャリア分布関数の値である分布関数値から、半導体デバイスの電気特性の基になるキャリア濃度分布を求めることができ、またこのキャリア濃度分布から、半導体デバイスの電圧−電流特性を始めとする、種々の電気特性を得ることができる。
【０００３】
ところで、このような方程式による演算処理にはコンピュータが利用されており、従来のシミュレーション装置にはポアソン方程式を解く演算処理部およびボルツマン方程式を解く演算処理部が設けられている。
両方程式は偏微分方程式で表わされ、特にボルツマン方程式は、時間、位置ベクトルおよび波数ベクトルの合計７つの変数の偏微分方程式で表わされることから、このボルツマン方程式を比較的容易に解くために、流体モデル法およびモンテカルロ法を利用することが提案されている（「プロセス・デバイスシミュレーション技術」、ｐ９１〜１３３およびｐ１９７〜２１４、檀良  編著、産業図書株式会社発行、昭和６３年）。
【０００４】
【発明が解決しようとする課題】
しかしながら、流体モデルによってボルツマン方程式を解く方法では、例えばキャリアの温度を仮定し、あるいはキャリアのエネルギー分布関数がマクスウェル−ボルツマンの分布関数である等、方程式の簡略化のために多数の仮定を導入する必要がある。そのため、この多数の仮定の導入により、現実の環境での結果と不一致を生じることがあり、常に導入した仮定についての妥当性が問題となる。
また、モンテカルロ法によってボルツマン方程式を解く方法では、方程式の簡略化およびそれに伴う多数の仮定の導入が行なわれることはないが、解法自体が乱数に基づくことから、確率の小さな現象の演算にも時間を要し、演算処理に多大の時間が掛かる。しかも、確率の小さな現象により、解そのものが引きずられ、このために現実には有り得ないような値を含む結果が求められることがあり、常に安定した結果を求めることが困難である。
【０００５】
そこで、半導体デバイスの電気特性を模擬的に、比較的速く、しかも現実に即する安定した結果として得ることのできる半導体デバイスシミュレーション方法および装置の現出が望まれていた。
【０００６】
【課題を解決するための手段】
本発明は、基本的には、ボルツマン方程式を解くにあたり、いわゆるセルオートマトン（ｃｅｌｌｕｌａｒ  ａｕｔｏｍａｔｏｎ，  自己増殖機械）法を利用する。
先ず、本発明の基本概念を述べる。
ボルツマン方程式では、分布関数が時間、位置ベクトルおよび波数ベクトルの関数として表わされる。
【０００７】
これに対応して、本発明では、位置ベクトルに係わる位置空間および波数ベクトルに係わる運動量空間をそれぞれ生成する。
位置空間を小領域に区画し、各小領域にそれぞれの小領域におけるキャリアの存在確率を示すキャリアの分布関数値を代表する代表点を設ける。また、位置空間の各代表点毎にそれぞれ運動量空間を対応させ、各運動量空間を小領域に区画し、運動量空間の各小領域にそれぞれの小領域におけるキャリアの存在確率を示すキャリアの分布関数値を代表する代表点を設ける。
各位置空間の代表点に初期値を与え、この初期値に基づき対応する運動量空間に初期値を与える。また、この位置空間の初期値に対応してキャリアを各運動量空間の各小領域に配分する。
【０００８】
各運動量空間でキャリアが散乱されるまでの時間Δｔを求め、後述するボルツマン方程式の衝突までの走行に関する項についての解を求めるために、運動量空間の時間Δｔ後の運動量空間での運動量の変化量を求め、さらに、この運動量空間に対応する位置空間での位置の変化量を求める。
これら両変化量から、位置空間および運動量空間における各小領域でのキャリア存在確率すなわち各代表点における分布関数値が補正され、時間Δｔ後の各小領域でのキャリア存在確率が求められる。
【０００９】
また、ボルツマン方程式の散乱に関する項についての解を求めるために、運動量空間においてキャリアの散乱後の運動量を求め、この運動量空間における散乱後の運動量の変化によるキャリア存在確率である分布関数値を求め、この運動量空間の分布関数値から対応する位置空間での分布関数値を求める。
さらに、この位置空間および運動量空間の各小領域でのキャリア存在確率を示す各分布関数値を初期値として、再び、散乱されるまでの時間Δｔを求め、時間Δｔ後の走行に関する項についての解および散乱に関する項についての解を求めるステップが所定回数繰り返される。
以上がボルツマン方程式を解くにあたっての本発明の方法および装置についての基本概念である。
【００１０】
本発明に係る半導体デバイスシミュレーション方法は、上記の基本概念に沿って半導体デバイスの電気的特性をコンピュータを用いて模擬的に求める半導体デバイスシミュレーション方法において、対象とする半導体デバイスの不純物濃度分布および境界条件についての条件を入力し、入力条件に沿って、電位とキャリア濃度との関係を示すポアソン方程式を解いて電界に関する情報を求め、ポアソン方程式を解くことにより求められた電界に基づき、ボルツマン方程式を解いてキャリアの分布関数を求め、ボルツマン方程式を解くことにより求められたキャリアの分布関数値から求められるキャリア濃度をポアソン方程式に代入して両方定式の演算を繰り返し、集束する分布関数に基づいて対象とする半導体デバイスの電気特性を求める。
【００１１】
ボルツマン方程式の解として、以下のステップ
１．材料の形状に対応してメッシュ状に区画された多数の小領域からなる各位置空間を生成し、該各位置空間に代表点を設け、該各代表点のそれぞれにキャリア存在数を示す初期値を与え、該初期値を前記分布関数Ｆに代入して分布関数値（Ｆ）を算出し、
２．前記位置空間の各代表点毎にメッシュ状に区画された多数の小領域からなる運動量空間を生成し、該各運動量空間に代表点を設け、該代表点のそれぞれにキャリア存在数を示す初期値を与え、該初期値を前記分布関数値（Ｆ）に代入して分布関数値（ｆ）を算出し、
３．前記運動量空間の分布関数値（ｆ）に対応してキャリアを対応する運動量空間に配分し、
４．半導体材料の物理モデルに応じて該運動量空間に配分されたキャリアが散乱されるまでの時間Δｔを求め、
５．ポアソン方程式の解から求められた電界及びエネルギーバンド構造に基づいて全運動量空間のキャリアの加速運動による時間Δｔ後の位置および運動量についての各変化量を求め、
６．キャリアの運動後の前記位置空間および前記運動量空間上の位置を基に、各空間の各小領域の前記代表点における分布関数値（Ｆ、ｆ）を補正し、
７．前記運動量空間において、散乱の種類に応じたキャリア散乱後の運動量を求め、当該運動量空間上の位置を基に、該運動量空間の各小領域の前記各代表点における分布関数値（ｆ）を補正し、
８．分布関数値（Ｆ、ｆ）が補正された前記位置空間および前記運動量空間について、前記３．から前記７．のステップを所定回数繰り返し、分布関数値（Ｆ）から求められてポアソン方程式に代入されるキャリア濃度を逐次更新する、
の実行により求められることを特徴とする。
【００１２】
また、本発明に係る半導体デバイスシミュレーション装置は、対象とする半導体デバイスの不純物濃度分布および境界条件を含む条件が入力される入力装置と、該入力装置からの入力条件に基づいてキャリア分布関数についての演算を行うプロセッサと、該プロセッサから出力されるキャリア分布関数に関する情報を表示する出力装置とを含む。
【００１３】
プロセッサは、入力装置からの条件に沿って電界に関するポアソン方程式を解く演算処理部と、ポアソン方程式を解くことにより求められた電界に基づき、キャリアの分布関数を求めるボルツマン方程式を解く演算処理部と、両演算処理部の反復演算処理により求められたキャリア分布関数が集束するか否かを判定する集束判定部とを備える。
【００１４】
ボルツマン方程式を解く演算処理部は、メッシュ状に区画された小領域からなる位置空間を生成する位置空間生成部分と、位置空間の各小領域に対応する複数の運動量空間であってそれぞれがメッシュ状に区画された小領域からなる運動量空間を生成する運動量空間生成部分と、位置空間の各小領域に与えられたキャリアの分布関数値（Ｆ）に基づく分布関数値（ｆ）が与えられた運動量空間に位置空間での分布関数値（Ｆ）に対応して配分されたキャリアが散乱されるまでの時間Δｔを求めるモデル演算処理部分と、該演算処理部分により求められた時間Δｔ後のキャリアの位置空間および運動量空間における位置および運動量についての変化量を求め、該変化量に基づいて位置空間および運動空間の各小領域の分布関数値（Ｆ、ｆ）を補正する走行項演算処理部分と、位置空間の各小領域の補正された分布関数値（Ｆ）に対応する運動量空間において、散乱の種類に応じたキャリア散乱後の運動量を求め、当該運動量空間上の位置を基に、該運動量空間の各小領域の分布関数値（ｆ）を補正する散乱項演算処理部分とを有する。
【００１５】
集束判定部は、走行項演算処理部分および散乱項演算処理部分により補正された位置空間および運動量空間における各小領域での分布関数値（Ｆ、ｆ）をもとに、モデル演算処理部分、走行項演算処理部分および散乱項演算処理部分の所定回数の反復演算処理およびポアソン方程式を解く演算処理部への反復演算処理を経て得られる分布関数値（Ｆ）から求められるキャリアの分布関数が集束するか否かを判定する。
【００１６】
【作用】
本発明の方法では、位置空間の各代表点に与えられる初期値は、従来の流体モデル法におけると同様に、仮定的に導入されるが、この流体モデル法におけるように仮定条件が固定的に取り扱われることはなく、ボルツマン方程式の反復演算およびポアソン方程式の演算をも含む反復演算によって、これら仮定的に導入される条件等は、適性値に向けて順次更新される。従って、従来のように特定の仮定についての妥当性が問題となることはない。
また、乱数に基づく演算処理ではないことから、演算処理にモンテカルロ法におけるような多大の演算処理時間を要することはなく、しかも、確率の小さな現象による解のばらつきが防止される。
【００１７】
また、本発明の装置では、該装置を構成する各部をコンピュータにより構成することにより、本発明の方法を適切かつ有効に実施することができる。これにより、従来に比較して高い精度でしかも高速に適正な半導体デバイス特性のシミュレーション結果を得ることができる。
【００１８】
【実施例】
以下、本発明を図示の実施例に沿って詳細に説明する。
図１は、本発明に係るシミュレーション装置のブロック図である。
本発明に係る半導体デバイスシミュレーション装置１０は、対象とする半導体デバイスの構造的、電気的諸条件および必要とされる初期条件等が入力される入力装置１１と、該入力装置からの情報に基づいて半導体デバイスのキャリア分布関数についての演算を行うプロセッサ１２と、該プロセッサから出力されるキャリア分布関数に関する情報を表示する出力装置１３とを含み、プロセッサ１２には、その演算処理のために必要な情報を格納するメモリ１４が接続されている。
【００１９】
入力装置１１から入力される構造的あるいは電気的条件としては、対象とする半導体デバイスが例えばＭＯＳトランジスタのような素子であるとき、ゲート長、ゲート幅、ゲート絶縁膜の厚さのような形状に関する条件、半導体基板にソースあるいはドレイン領域を形成するためにこの基板に導入される不純物の濃度分布、接合部における境界条件等の解析に必要な条件がある。また初期条件として、キャリア濃度（ｐ−ｎ：正孔と電子の差）が不純物濃度（ＮＤ  −ＮＡ  ：ドナー濃度とアクセプタ濃度との差）に等しいと仮定する電荷中性条件のような初期値条件が入力される。
【００２０】
プロセッサ１２は、入力装置１１からの条件に沿って、電位とキャリア濃度の関係を示すポアソン方程式（１式）を解く演算処理部１５と、ポアソン方程式を解くことにより求められた電界に基づき、キャリアの分布関数を求めるボルツマン方程式（２式）を解く演算処理部１６と、両演算処理部１５、１６の反復演算処理により求められたキャリア分布関数が集束するか否かを判定する集束判定部１７とを備える。
【００２１】
ポアソン方程式の演算処理部１５は、例えば初期条件として入力された前記電荷中性条件に沿った電位を算出する。この電位は、引き続くボルツマン方程式の演算処理部１６のために微分により、電界に変換される。
【００２２】
この電界に基づいてキャリア分布関数を求める演算処理部１６は、メッシュ状に区画された小領域からなる位置空間を生成する位置空間生成部分１８と、位置空間の各小領域に対応する複数の運動量空間であってそれぞれがメッシュ状に区画された小領域からなる運動量空間を生成する運動量空間生成部分１９と、運動量空間でキャリアが散乱されるまでの時間Δｔを求めるモデル演算処理部分２０と、ボルツマン方程式（２式）の左辺で示されるキャリア走行項についての演算処理を行う走行項演算処理部分２１と、ボルツマン方程式（２式）の右辺で示されるキャリア散乱項についてキャリアの散乱後の運動量を求める散乱項演算処理部分２２とを有する。
【００２３】
ボルツマン方程式（２式）に示されるキャリア分布関数Ｆは、位置についての位置ベクトルおよび運動量についての波数ベクトルの関数であり、これに関連して、位置空間生成部分１８および運動量空間生成部分１９は、それぞれ位置空間および運動量空間をメモリ１４に生成する。
【００２４】
図２は位置空間を概念的に示す説明図であり、図３および図４はそれぞれ運動量空間を概念的に示す説明図である。
図示の例では、位置空間２３および運動量空間２４は、説明の簡略化のために、それぞれ直交するメッシュによって多数の小領域に区画された２次元メッシュで表現されている。これら空間２３および２４を３次元で表わすことができる。
また、位置空間２３の実線は、キャリア分布関数を求めようとする材料の形状に対応している。
【００２５】
位置空間２３の各小領域は、キャリアの存在確率である分布関数値（Ｆ）をそれぞれ与えられ、各小領域の分布関数値（Ｆ）を代表する代表点が例えばメッシュの交点あるいは小領域の中央部等にそれぞれ設定される。図２には、図面の簡略化のために小領域２３Ａの分布関数値（Ｆ）を代表する１つのキャリアＣ１のみが示されている。
位置空間２３の１つのキャリアＣ１に対して、例えば図３で示したような１つの運動量空間２４が対応する。従って、位置空間２３の小領域の数に対応した数の運動量空間２４が生成される。
【００２６】
各運動空間２４の小領域には、該運動空間に対応する位置空間２３のキャリアＣ１の分布関数値（Ｆ）に対応した分布関数値（ｆ）が与えられる。例えば、キャリアＣ１の分布関数値（Ｆ）が１であったとすると、このキャリアＣ１に対応する運動量空間２４の全小領域における分布関数値（ｆ）の和すなわち対応する運動量空間２４の全キャリアＣ２の分布関数値（ｆ）の和が、対応するキャリアＣ１の分布関数値（Ｆ）である１に等しくなるように設定される。
【００２７】
この運動量空間２４のキャリアＣ２の数については、複数のキャリアＣ２にそれぞれ等価の重みを与えるか、あるいは単一のキャリアＣ２にこれに対応するキャリアＣ１の分布関数値（Ｆ）に応じた分布関数値（ｆ）を与えることができる。
複数のキャリアＣ２に等価の重みを与える場合、例えばキャリアＣ１の分布関数値（Ｆ）が１であり、これに対応するキャリアＣ２が図３に示されているように２個であるとき、各キャリアＣ２には、それぞれ０．５の重みが与えられる。従って、０．５の重みが与えられたキャリアＣ２は、これに与えられた分布関数値（ｆ）の半値の存在確率を持つことになる。
【００２８】
図４には運動量空間２４の単一のキャリアＣ２が示されており、この単一のキャリアＣ２には１．０の重みが与えられている。従って、この単一のキャリアＣ２は、対応するキャリアＣ１の分布関数値（Ｆ）に応じて与えられた分布関数値（ｆ）を持つこととなる。
次に述べるキャリア走行項についての運動量空間２４で各キャリアの運動を追う必要があることから、このように各運動量空間２４毎に単一のキャリアＣ２を適用することが、精度の低下を招くことなく演算処理時間の短縮化を図る上で、有利である。
【００２９】
ボルツマン方程式（２式）の左辺で示されるキャリア走行項は、キャリアが衝突すなわち散乱するまでの時間Δｔの間におけるキャリアの走行運動を表わす項である。この時間Δｔを求めるモデル演算処理部分２０は、半導体材料に応じたキャリアの物理モデルに応じて、例えば格子散乱、不純物散乱等毎に散乱が生じるまでの時間Δｔを算出する。
【００３０】
また、ボルツマン方程式（２式）の右辺である散乱項を零におき、この式に、ポアソン方程式演算処理部１５で求められた電界およびモデル演算処理部分２０で求められた時間Δｔを代入して演算処理することにより、散乱を生じるまでのキャリアの運動量空間２４での運動を求めることができる。
走行項演算処理部分２１は、キャリアに応じたエネルギバンドモデルを用いて、運動量空間２４でのキャリアＣ２の運動量の変化量を求め、また、求めた運動量空間２４での変化量をもとに位置空間２３での対応するキャリアの位置についての変化量を求める。
このキャリアの位置の変化量は、エネルギーバンドモデルと電界とから求めることもできる。
【００３１】
図４には、運動量空間２４で黒丸位置から白丸位置に移動したキャリアＣ２が示されている。走行項演算処理部分２１は、このキャリアＣ２の変化量を算出すると共に、この移動によって生じる運動量空間２４の各小領域での分布関数値（ｆ）を補正する。
【００３２】
この補正については、図４を参照するに、キャリアＣ２の分布関数値（ｆ）は、移動後のキャリアＣ２を取り巻く４つの代表点Ａ、Ｂ、Ｃ、ＤおよびキャリアＣ２の移動後の点で区画される面積Ｓ１、Ｓ２、Ｓ３、Ｓ４に応じて、次式
Ａ＝Ｓ３／Ｓ
Ｂ＝Ｓ４／Ｓ
Ｃ＝Ｓ１／Ｓ
Ｄ＝Ｓ２／Ｓ
（各式において、ＳはＳ１，Ｓ２，Ｓ３，Ｓ４の和である。）
で示される割合で、近接する各代表点Ａ、Ｂ、Ｃ、Ｄに分け与えられる。
従って、走行項演算処理部分２１は、各代表点をそれぞれの新たな分布関数値（ｆ）に補正する。
なお、分布関数値（ｆ）は、図４に沿って説明したこの補正に先立って、これと同様の方法により、キャリアの位置空間での移動に応じて、走行項演算処理部分２１により予め補正されている。
【００３３】
位置空間２３あるいは運動量空間２４が１次元で表わされている場合、分布関数値（Ｆ、ｆ）の補正は、先に述べた面積比におけると同様に、近接する代表点間での距離比に応じて補正され、また位置空間２３あるいは運動量空間２４が３次元で表わされている場合、分布関数値（Ｆ、ｆ）の補正は、体積比に応じて同様に補正される。
また、図３に示したように、複数の等価のキャリアＣ２が適用されている場合、走行項演算処理部分２１は、各キャリアＣ２が移動先で取り囲まれた、それぞれ斜線で示す小領域の代表点に、分布関数値（ｆ）をキャリアの重みで等分した値を与える。
【００３４】
ボルツマン方程式（２式）の右辺で示されるキャリア散乱項についてキャリアの散乱後の運動量を求める散乱項演算処理部分２２は、運動量空間２４で生じる各種の散乱についてその散乱後のキャリアの運動量を求め、この運動量空間２４での運動量に基づいて、走行項演算処理部分２１におけると同様に、運動量空間２４の各小領域での代表点の分布関数値（ｆ）を補正する。
キャリアの散乱後の運動量を求めるについて、散乱の種類を確率分布に基づく乱数に応じて選択することができる。
しかしながら、より迅速かつ高精度の演算処理を求める上で、各散乱モデルについての確率を求め、この確率に基づいてキャリア運動量を算出することが望ましい。
【００３５】
図５は、運動量空間２４におけるキャリアＣ２の散乱での運動量変化による分布関数値（ｆ）の補正についての説明図である。
黒丸で示されるキャリアＣ２が確率Ｐ１で生じる散乱により白丸の位置（Ｃ２−１）に移動し、また確率Ｐ２で生じる他の散乱により２重の白丸の位置（Ｃ２−２）に移動したとする。このとき、キャリアＣ２−１には、キャリアＣ２が有していた分布関数値（ｆ）に確率Ｐ１の重みを掛けた値が与えられ、またキャリアＣ２−２には、キャリアＣ２が有していた分布関数値（ｆ）に確率Ｐ２の重みを掛けた値が与えられる。
【００３６】
このキャリアの散乱による補正は、図４について説明したと同様に、キャリアＣ２−１については、キャリアＣ２の分布関数値（ｆ）に確率Ｐ１を掛けた値が、キャリアＣ２−１を取り巻く４つの代表点Ａ、Ｂ、Ｅ、Ｆにそれぞれ区画面積Ｓ１、Ｓ２、Ｓ３、Ｓ４に応じて、次式
Ａ＝Ｓ４／Ｓ
Ｂ＝Ｓ３／Ｓ
Ｅ＝Ｓ２／Ｓ
Ｆ＝Ｓ１／Ｓ
（各式において、ＳはＳ１，Ｓ２，Ｓ３，Ｓ４の和である。）
で示される割合で、近接する各代表点Ａ、Ｂ、Ｅ、Ｆに分け与えられる。
【００３７】
また、キャリアＣ２−２についても、キャリアＣ２の分布関数値（ｆ）に確率Ｐ２を掛けた値が、キャリアＣ２−２を取り巻く４つの代表点Ｂ、Ｃ、Ｄ、Ｅに同様な区画面積比（Ｓ１，Ｓ２，Ｓ３，Ｓ４）に応じて分け与えられる。
これらの分布関数値（ｆ）の配分により、運動量空間２４の分布関数値（ｆ）の散乱による補正がなされ、この運動量空間２４における各代表点の分布関数値の補正により、ボルツマン方程式の基本的な１サイクルの演算処理が終る。
【００３８】
この１サイクルの演算処理後、走行項演算処理部分２１により求められた位置空間２３における各分布関数値および散乱項演算処理部分２２により求められた運動量空間２４における各分布関数値を、それぞれ両空間２３および２４の初期値として、モデル演算処理部分２０、走行項演算処理部分２１および散乱項演算処理部分２２での前記したと同様な演算処理サイクルが所定回数繰り返される。
【００３９】
モデル演算処理部分２０、走行項演算処理部分２１および散乱項演算処理部分２２の反復演算が所定回数に達すると、集束判定部１７が走行項演算処理部分２１により求められた位置空間２３における各分布関数値（Ｆ）から求められるキャリアの分布関数Ｆが集束したか否か、すなわち分布関数Ｆがほぼ誤差内で集束したか否かを判定する。
集束していなければ、その分布関数値（Ｆ）で求められるキャリア濃度を演算処理部１５におけるキャリア濃度として、再度、各演算処理部１５、１６での演算処理が繰り返される。
【００４０】
出力装置１３は、この集束するキャリア分布関数Ｆを解として表示する。この出力装置１３に、キャリア分布関数Ｆから求められる種々の電気特性に関する情報を必要に応じて表示させることができる。
【００４１】
図６および図７は、本発明のシミュレーション方法の手順を示すフローチャートである。
半導体デバイスシミュレーション装置１０の作動の概略を図６に沿って説明する。
【００４２】
先ず、入力装置１１に、対象とする半導体デバイスの不純物濃度分布、境界条件あるいは電荷中性条件等の初期条件が設定値設定として入力される（ステップＳ６１）。
これらの入力条件に沿って、ポアソン方程式を解く演算処理部１５が演算処理を行い、電界に関する情報を求める（ステップＳ６２）。
続いて、ポアソン方程式を解く演算処理部１５により求められた電界に基づき、ボルツマン方程式を解く演算処理部１６がセルオートマトン法に沿って、キャリア分布関数Ｆを求める（ステップＳ６３）。
演算処理部１６により求められたキャリア分布関数Ｆが、集束するか否かを集束判定部１７により、判定される（ステップＳ６４）。
【００４３】
キャリア分布関数Ｆが集束していない場合、このキャリア分布関数Ｆから求められるキャリア濃度が演算処理部１５のキャリア濃度として、再度ステップＳ６２以下の演算処理が繰り返される。
また、キャリア分布関数Ｆが集束する場合、出力装置１３により、キャリア分布関数Ｆに関する情報が出力される。
【００４４】
先に述べたボルツマン方程式を解く演算処理（ステップＳ６３）では、図７に示されているように、位置空間生成部分１８および運動量空間生成部分１９によって位置空間および運動空間がそれぞれ生成されかつそれぞれがメッシュ状に区画され、各空間の小領域について代表点が決められ、位置空間の個々の代表点についての分布関数値（Ｆ）が初期値として設定され、これに基づいて運動量空間の各初期値が決まる（ステップＳ７１）。
【００４５】
この位置空間の初期値に、マクスウェル−ボルツマン分布関数を適用することができ、また、従来の前記した流体モデル法によって得られた分布関数を適用することができる。この流体モデル法によって得られた分布関数は、適正な分布関数に近似する関数であることから、これに基づき各代表点の初期値を設定することにより、すなわち、セルオートマトン法の初期解に流体モデル法を用いることにより、演算処理時間の短縮化を図ることができる。
【００４６】
ステップＳ７１に示される初期設定は、図７に示されたサブルーチンを実行する際の最初に実行されるのみであり、集束ステップＳ６４からの繰り返しによっる二回目以降のループではこのステップＳ７１が繰り返されることはない。
【００４７】
また、運動量空間生成部分１９は、運動量空間の分布関数値（ｆ）に対応して、キャリアを対応する運動量空間２４に配分する（ステップＳ７２）。
運動量空間２４に配分されたキャリアが散乱するまでの時間Δｔがモデル演算処理部分２０により求められ（ステップＳ７３）、この時間Δｔ後のキャリアの位置空間２３および運動量空間２４での変化量が走行項演算処理部分２１により求められ（ステップＳ７４）、その結果に基づいて位置空間２３および運動量空間２４での分布関数値（Ｆ、ｆ）が補正される（ステップＳ７５）。
【００４８】
続いて、散乱項演算処理部分２２が運動量空間２４でのキャリア散乱後の運動量を求め、この運動量に基づいて、運動量空間２４の分布関数値（ｆ）を補正する（ステップＳ７６）。
この補正された位置空間２３および運動量空間２４での分布関数値（Ｆ、ｆ）に基づいて、ステップＳ７２〜ステップＳ７６が所定回数反復され（ステップＳ７７）、その後、図６に示した集束判断ステップＳ６４に移る。
【００４９】
本発明の半導体デバイスシミュレーション装置１０によれば、位置空間の各代表点に与えられる初期値を決める仮定条件は、固定的に取り扱われることはなく、ボルツマン方程式の反復演算およびポアソン方程式の演算をも含む反復演算によって、適性値に向けて順次更新される。従って、従来のように特定の仮定についての妥当性が問題となることはない。
また、乱数に基づく演算処理ではないことから、演算処理にモンテカルロ法におけるような多大の演算処理時間を要することはなく、しかも、確率の小さな現象による解のばらつきが防止される。
従って、迅速に、しかも高い精度で適正なシミュレーション結果を得ることができる。
【００５０】
図８は、半導体デバイスシミュレーション装置１０によって得られたキャリア濃度分布特性曲線を示すグラフである。
この例では、基板２５としてＰ型シリコンを用いたＭＯＳＦＥＴ２６について、ゲート電極２７の両側にソース／ドレイン領域２８を形成するために、例えばＡｓ＋のようなイオンが注入され、ある外部電圧が電極に印加されたときのキャリア濃度の分布特性曲線がＭＯＳＦＥＴ２６の縦断面図上に求められた。
それぞれ１０１８〜１０２０／ｃｍ３  の等濃度を示す曲線が滑らかな特性曲線で表わされているように、モンテカルロ法におけるような統計的揺らぎによる特異な点が出現することなく、より正確なキャリア濃度分布特性曲線が迅速に得られた。
【００５１】
図９は、運動量空間２４の他の例を示す説明図である。
運動量空間２４を、図３〜図５に示したように方形に区画することに代えて、図９に示すように、円形メッシュを適用することができる。図では説明の簡略化のために２次元で示されているが、デカルト座標では球形メッシュが適用される。
【００５２】
図９で、キャリアが散乱前の運動量を示す位置Ｃ３からある散乱によってＣ３−１に示した位置へ移動したとすると、この移動による分布関数値（ｆ）の補正は次のようになされる。
移動後の位置Ｃ３−１を取り囲む代表点Ａ、Ｂ、Ｃ、Ｄに分布関数値（ｆ）が配分されるが、先ずＸ軸およびＹ軸上への配分比が運動ベクトルのＸ軸成分およびＹ軸成分比であるｃｏｓθ、ｓｉｎθに沿って配分される。
【００５３】
次にＸ軸上に配分された分布関数値（ｆ）×ｃｏｓθの値が、Ｘ軸上の代表点Ｃ、Ｄ上に、図４に沿って説明したと同様な距離比に応じて、
Ｃ＝ｎ／（ｎ＋ｍ）
Ｄ＝ｍ／（ｎ＋ｍ）
の割合で配分される。
同様に、Ｙ軸上に配分された分布関数値（ｆ）×ｓｉｎθの値が、Ｙ軸上の代表点Ａ、Ｂ上に、その距離比に応じて代表点Ｃ、Ｄにおけると同様に配分される。
このＸ軸、Ｙ軸上の代表点の分布関数値（ｆ）の配分により、運動量空間２４の分布関数値（ｆ）の散乱による補正がなされる。
【００５４】
このように、円形メッシュあるいは球状メッシュを使用することにより、代表点を軸（Ｘ、Ｙ、またはＺ）上に集約することができ、演算の簡素化およびそれによる高速化を図ることができる。
【００５５】
また、以上に説明したところでは、各位置空間および運動量空間のメッシュで区画された各小領域はその領域内では均等な分布関数値すなわち存在確率を示すと考えたが、このような仮定で、演算時間の短縮化を図るためにメッシュを粗く設定すると、疑似拡散と呼ばれる現象が生じ、誤差が大きくなる。
これを防止するために、各小領域の存在確率を示す分布関数値を補間関数により修正することにより、メッシュを粗く設定しても疑似拡散を防止し、精度の高い演算をより短時間で行うことが可能になる。
【００５６】
図１０は、この補間関数の原理の説明図である。
横軸には、各小領域を区画するメッシュのＸ軸が示され、縦軸には各小領域の代表点が表わす存在確率が示されている。
各小領域の代表点がそれぞれ黒点で示されているが、これら各小領域が均等な存在確率を示すと仮定すると、分布関数は、図中破線で示すような折れ線で示されることとなる。しなしながら、このような仮定は非現実的であり、真の分布関数の近似に過ぎない。
そこで、各代表点を通る関数曲線を補間関数Ｈとして採用し、キャリアの分布関数値（Ｆ、ｆ）を隣接する各代表点に配分する際、この補間関数Ｈにより配分比を修正することにより、メッシュを比較的細かくすることなく疑似拡散現象を防止することができることから、より短時間での高精度の演算処置が可能になる。
【００５７】
本発明は以上の実施例に限定されない。位置空間および運動量空間の各Ｘ軸およびＹ軸をそれぞれ等間隔の等値で区画した例を示したが、各軸の間隔を不等とすることができる。この場合、円形メッシュは長円メッシュになる。
また、磁場を考慮することができ、この磁場についての演算処理過程をポアソン方程式の演算処理と、ボルツマン方程式の演算処理との間に、挿入することができる。また、磁場を考慮した場合には、走行項演算処理部分２１での演算に磁場の効果を考慮する必要がある。
【００５８】
【発明の効果】
以上説明した本発明の方法では、ボルツマン方程式を解くに際し、位置空間の各代表点に与えられる初期値等は仮定的に導入されるが、ボルツマン方程式の反復演算およびポアソン方程式の演算をも含む反復演算によって、これら仮定的に導入される諸条件は、固定的に取り扱われることはなく適性値に向けて順次更新される。そのため、従来のような特定の仮定についての妥当性が問題となることはない。また、多大の演算処理時間を要することはなく、確率の小さな現象による解のばらつきが防止される。
【００５９】
従って、本発明の方法によれば、現実に即する安定した適正なシミュレーション結果が迅速に得られる。
また、本発明の装置よれば、該装置を構成する各部をコンピュータにより構成することにより、本発明の方法を適切かつ有効に実施することができることから、適正なシミュレーション結果を、高い精度で、しかも迅速に得ることができる。
【図面の簡単な説明】
【図１】本発明に係るシミュレーション装置のブロック図である。
【図２】本発明に係る位置空間の説明図である。
【図３】本発明に係る運動量空間の説明図（その１）である。
【図４】本発明に係る運動量空間の説明図（その２）である。
【図５】本発明に係る運動量空間の説明図（その３）である。
【図６】本発明に係るシミュレーション方法の手順を示すフローチャート（その１）である。
【図７】本発明に係るシミュレーション方法の手順を示すフローチャート（その２）である。
【図８】本発明の装置によるシミュレーション結果を示すキャリア濃度分布特性曲線図である。
【図９】本発明に係る運動量空間の説明図（その４）である。
【図１０】本発明に係る補間関数の説明図である。
【符号の説明】
１０  半導体デバイスシミュレーション装置
１１  入力装置
１２  プロセッサ
１３  出力装置
１５  ポアソン方程式演算処理部
１６  ボルツマン方程式演算処理部
１７  集束判定部
１８  位置空間生成部分
１９  運動量空間生成部分
２０  モデル演算処理部分
２１  走行項演算処理部分
２２  散乱項演算処理部分
２３  位置空間
２４  運動量空間
【図１】
0003569750.tif 000002
【図２】
0003569750.tif 000003
【図３】
0003569750.tif 000004
【図４】
0003569750.tif 000005
【図５】
0003569750.tif 000006
【図６】
0003569750.tif 000007
【図７】
0003569750.tif 000008
【図８】
0003569750.tif 000009
【図９】
0003569750.tif 000010
【図１０】
0003569750.tif 000011