電子輸送のためのモンテカルロ法

ボルツマン輸送方程式モデルは、半導体における輸送現象の解析において主要なツールとして用いられてきた。BTE方程式は次のように表される。

$${\displaystyle {\frac {\partial f}{\partial t}}+{\frac {1}{\hbar }}\nabla _{k}E(k)\nabla _{r}f+{\frac {qF(r)}{\hbar }}\nabla _{k}f=\left[{\frac {\partial f}{\partial t}}\right]_{\mathrm {collision} }}$$

${\displaystyle v={\frac {1}{\hbar }}\nabla _{k}E(k)}$

分布関数f は、関心のあるすべての観測量を抽出するために使用される無次元関数であり、実空間とk 空間の両方で電子分布を完全に描写します。さらに、物理的には、位置rと時刻tにおけるエネルギーkの粒子の占有確率を表します。加えて、7 次元積分微分方程式 (位相空間に 6 次元、時間に 1 次元) であるため、BTE の解は複雑で、非常に特殊な制約の下で閉じた解析形式で解くことができます。数値的には、BTE の解は決定論的方法または確率的方法のいずれかを使用して使用されます。決定論的方法の解は、球面調和関数アプローチなどのグリッドベースの数値方法に基づいていますが、モンテカルロ法は BTE を解くために使用される確率的アプローチです。 

半古典モンテカルロ法は、複雑なバンド構造と散乱過程を含むボルツマン輸送方程式の厳密解を求めるために用いられる統計的手法です。この手法が半古典的である理由は、散乱機構はフェルミの黄金律を用いて量子力学的に扱われるのに対し、散乱事象間の輸送は古典的な粒子概念を用いて扱われるからです。モンテカルロモデルは、本質的に各自由飛行における粒子の軌跡を追跡し、対応する散乱機構を確率的に選択します。半古典モンテカルロ法の大きな利点の2つは、散乱項内の様々な異なる散乱機構を正確に量子力学的に扱うことができることと、エネルギー空間またはk空間におけるキャリア分布の形式に関する仮定がないことです。電子の運動を記述する半古典方程式は

${\displaystyle {\frac {dr}{dt}}={\frac {1}{\hbar }}\nabla _{k}E(k)}$

${\displaystyle {\frac {dk}{dt}}={\frac {qF(r)}{\hbar }}}$

ここで、F は電場、E(k) はエネルギー分散関係、k は運動量波数ベクトルです。上記の式を解くには、バンド構造 (E(k)) に関する深い知識が必要です。E(k) 関係は、状態密度( DOS) や粒子速度など、輸送に必要な有用な情報を示すだけでなく、デバイス内での粒子の動き方を記述します。半経験的擬ポテンシャル法を用いることで、フルバンド E(K) 関係を得ることができます。^{[ 1 ]}

ドリフト拡散（DD）モデルと流体力学（HD）モデルはどちらも、長チャネルデバイスに有効な簡略化された近似を使用して、ボルツマン輸送方程式（BTE）のモーメントから導出できます。DDスキームは最も古典的なアプローチであり、通常、ドリフトと拡散の成分を考慮して、キャリアのポアソン方程式と連続方程式を解きます。このアプローチでは、電荷の通過時間はエネルギー緩和時間に比べて非常に大きいと仮定されます。^{[ 2 ]}一方、HD法は、BTEのモーメントから得られるエネルギーバランス方程式を使用してDDスキームを解きます。^{[ 3 ] [ 4 ]}これにより、キャリアの加熱や速度オーバーシュート効果などの物理的な詳細を捉えて計算することができます。言うまでもなく、支配方程式は強く結合しており、DDスキームと比較してより多くの変数を扱う必要があるため、HDシミュレーションでは正確な離散化方法が必要です。

半古典モデルの精度は、トランジスタ構造における重要な短チャネル効果(SCE) である古典的な速度オーバーシュート問題をどのように扱うかを調べることにより、BTE に基づいて比較されます。本質的に、速度オーバーシュートは、電流駆動と相互コンダクタンスの実験的に観測される増加に関連する、スケーリングされたデバイスの非局所効果です。 ^{[ 5 ]}チャネル長が小さくなると、速度は高電界領域で飽和しなくなり、予測される飽和速度をオーバーシュートします。この現象の原因は、キャリアの通過時間がエネルギー緩和時間と同程度になり、そのため、移動キャリアが短チャネルデバイスでの散乱によって印加電界と平衡に達するのに十分な時間がないことです。^{[ 6 ]} DD モデルと HD モデルを使用したシミュレーション結果 (Illinois Tool: MOCA) の概要を隣の図に示します。図 (a) では、チャネル領域全体で速度オーバーシュート効果を引き起こすほど電界が高くない場合を示しています。この限界では、DD モデルのデータはオーバーシュートのない領域で MC モデルによく適合しますが、HD モデルはその領域で速度を過大評価します。速度オーバーシュートは MC データではドレイン接合付近でのみ観測され、HD モデルはその領域でよく適合します。MC データから、速度オーバーシュート効果は高電界領域で急激であることがわかりますが、これは HD モデルには適切に含まれていません。図 (b) に示す高電界条件では、速度オーバーシュート効果はチャネルのほぼ全体に及び、チャネル領域では HD 結果と MC 結果が非常に近い値になります。

バンド構造は、エネルギー(E)と波数ベクトル(k)の関係を表します。バンド構造は、電場の作用下におけるキャリアの動き、散乱率、および衝突後の最終状態を計算するために使用されます。シリコンのバンド構造とそのブリルアンゾーンは下図に示されていますが、ブリルアンゾーン全体を満たす解析式は存在しません。いくつかの近似を用いることで、バンド構造には放物線モードと非放物線モードという2つの解析モデルが存在します。

電子がt付近のdtの間に次の衝突を起こす確率は、次式で与えられる。

${\displaystyle p(t)\,dt=P[k(t)]\exp[-\int _{0}^{t}P[k(t')]\,dt']\,dt}$

ここで、P[k(t)]dt は、状態 k の電子が時間 dt の間に衝突を起こす確率です。指数部の積分の複雑さのため、上記の式の分布で確率的自由飛行を生成することは現実的ではありません。この困難を克服するために、人々は仮想的な「自己散乱」スキームを使用します。これにより、この自己散乱を含む全散乱率は一定となり、例えば、${\displaystyle \Gamma }$ランダム選択により、自己散乱が選択された場合、衝突後の k′ は k と同じになり、キャリアは摂動を受けることなく飛行を続ける。定数を導入する${\displaystyle P(k)=\tau _{0}^{-1}}$上記の式は、

${\displaystyle p(t)={\frac {1}{\tau _{0}}}\exp(-t/\tau _{0}).}$

乱数rは非常に簡単に使用して確率的な自由飛行を生成することができ、その飛行時間は次のように与えられます。${\displaystyle t_{r}=-\tau _{0}\ln(r)}$自己散乱に使用されるコンピュータ時間は、自由飛行時間の計算の簡略化によって十分に補われます。^{[ 9 ]}自由飛行時間の計算速度を向上させるために、「定数法」や「区分的法」などのいくつかの手法が自己散乱イベントを最小限に抑えるために使用されています。

半導体デバイスの重要な電荷輸送特性、例えばオームの法則からの逸脱やキャリア移動度の飽和などは、散乱機構の直接的な結果である。したがって、半導体デバイスのシミュレーションにおいて、このような機構の物理を捉えることは非常に重要である。この点において、半導体モンテカルロシミュレーションは、ほぼ網羅的な散乱機構を容易かつ高精度に組み込むことができるため、非常に強力なツールである。自由飛行の持続時間は散乱率から決定される。各飛行の終わりに、散乱キャリアの最終エネルギー、あるいはそれと同等に、その新しい運動量と散乱角を決定するために、適切な散乱機構を選択する必要がある。この意味で、2つの物体間の衝突に関する古典的な運動論から自然に導き出される2つの大きなタイプの散乱機構を区別することができる。

弾性散乱とは、粒子が散乱された後もエネルギーが保存される現象である。したがって、弾性散乱では粒子の運動量の方向のみが変化する。不純物散乱と表面散乱は、おおまかに言えば、弾性散乱過程の良い例である。

非弾性散乱では、散乱粒子と散乱中心の間でエネルギーが伝達されます。電子フォノン相互作用は、散乱粒子によって特定のエネルギーのフォノンが放出または吸収されるため、本質的に非弾性です。散乱メカニズムをより数学的に詳細に特徴付ける前に、半導体モンテカルロシミュレーションを実行する際には、主に次のタイプの散乱イベントを扱う必要があることに注意することが重要です。^{[ 9 ]}

音響フォノン：電荷キャリアは、結晶格子内の原子の振動の音響モードとエネルギーを交換する。音響フォノンは主に結晶格子の熱励起によって発生する。

極性光学：電荷キャリアは結晶格子の極性光学モードのいずれかとエネルギーを交換します。これらのモードは共有結合半導体には存在しません。光学フォノンは、最小単位セルに複数の原子が存在する場合に、異なる種類の原子同士の振動によって発生し、通常は光によって励起されます。

非極性光学フォノン：エネルギーは光学モードと交換される。非極性光学フォノンは、一般的に共有結合性半導体およびGaAsのLバレーにおいて考慮する必要がある。

等価インターバレーフォノン：フォノンとの相互作用により、電荷キャリアは初期状態から、異なるが等価なバレーに属する最終状態に遷移します。通常、このタイプの散乱メカニズムは、電子が1つのXバレーから別のXバレーへ、または1つのLバレーから別のLバレーへ遷移することを記述します。^{[ 10 ]}

非等価インターバレーフォノン：異なる種類のバレー間での電荷キャリアの遷移を伴う。

圧電フォノン：低温用。

イオン化不純物:結晶格子中のイオン化不純物とのクーロン相互作用により、粒子が弾道軌道から逸脱することを反映しています。電子の質量は不純物の質量に比べて比較的小さいため、クーロン断面積は初期状態と最終状態の運動量の差とともに急速に減少します。^{[ 9 ]}したがって、不純物散乱事象は主に谷内散乱、バンド内散乱、そしてわずかな範囲でバンド間散乱について考慮されます。

キャリア間: (電子-電子、正孔-正孔、電子-正孔の相互作用)。キャリア濃度が高い場合、このタイプの散乱は電荷キャリア間の静電相互作用を反映します。この問題は、アンサンブルシミュレーション内の粒子数の増加に伴い、非常に計算負荷が高くなります。この範囲では、粒子とその周囲の電荷ガスとの短距離および長距離の相互作用を区別する粒子-粒子-粒子メッシュ (P3M) アルゴリズムが、半導体モンテカルロシミュレーションにキャリア間相互作用を含めるのに効果的であることが証明されています。^{[ 11 ]}多くの場合、キャリアの電荷はクラウドインセル法を使用してグリッドに割り当てられます。この方法では、特定の粒子の電荷の一部が、特定の重み係数で、最も近いグリッドポイントの特定の数に割り当てられます。

プラズモン：電荷キャリアの集団振動が特定の粒子に及ぼす影響を反映する。

モンテカルロシミュレーションに散乱を組み込むための計算効率の良いアプローチは、個々の機構の散乱率をテーブルに格納することである。特定の粒子状態に対する異なる散乱率が与えられれば、自由飛行の終わりに散乱過程をランダムに選択することができる。これらの散乱率は、多くの場合、ボルン近似を用いて導出される。ボルン近似では、散乱事象は、関与するキャリアの2つの運動量状態間の遷移にすぎない。セクションII-Iで議論したように、キャリアとその周囲環境（フォノン、電子、正孔、プラズモン、不純物など）との相互作用から生じる量子多体問題は、準粒子近似を用いて2体問題に還元できる。準粒子近似では、対象となるキャリアを結晶の残りの部分から分離する。^{[ 9 ]}これらの近似の範囲内では、 フェルミの黄金律は、一次近似で、ある状態からの散乱機構の単位時間あたりの遷移確率を与える。${\displaystyle |k\rangle }$州へ${\displaystyle |k'\rangle }$:

${\displaystyle S(k,k')={\frac {2\pi }{\hbar }}\left|\langle k|H'|k'\rangle \right|^{2}\cdot \delta (E-E')}$

ここで、H' は衝突を表す摂動ハミルトニアンであり、E と E′ はそれぞれ、キャリアと電子およびフォノンガスの両方から構成されるシステムの初期エネルギーと最終エネルギーである。ディラック${\displaystyle \delta }$-関数はエネルギー保存を表します。さらに、この用語は${\displaystyle \langle k|H'|k'\rangle }$一般に行列要素と呼ばれるものは、搬送波の初期波動関数と最終波動関数の内積を数学的に表す。^{[ 12 ]}

${\displaystyle \langle k|H'|k'\rangle ={\frac {1}{Vol}}\int _{\mathrm {Vol} }\psi _{k}(r)H'\psi _{k'}^{*}(r)\,dr}$

結晶格子では、波動関数は ${\displaystyle \psi _{k}(r)}$そして ${\displaystyle \psi _{k'}(r)}$これらは単にブロッホ波である。可能な場合、不純物散乱^{[ 13 ]}や音響フォノン散乱^{[ 14 ]}の場合のように、ハミルトニアンH'をフーリエ展開することで、行列要素の解析的表現が一般的に求められる。波数ベクトルqと周波数を持つフォノンによるエネルギー状態Eからエネルギー状態E'への遷移という重要なケースでは、${\displaystyle \omega _{q}}$エネルギーと運動量の変化は次のようになります。

${\displaystyle E'-E=E(k')-E(k)\pm \hbar \omega _{q}\,}$

${\displaystyle k'-k\pm q={\begin{cases}0&{\text{ }}\\R&{\text{Umklapp-process}}\end{cases}}}$

ここで、Rは逆格子ベクトルです。ウムクラップ過程（または U 過程）は、散乱後に粒子の運動量を変化させるため、半導体結晶の伝導を制限します。物理的には、U 過程は、粒子の最終運動量が第一ブリルアンゾーンの外側を向いているときに発生します。状態 k から状態 k' への単位時間あたりの散乱確率がわかれば、特定の散乱過程の散乱率を決定することは興味深いことです。散乱率は、逆格子空間内の状態kから他の任意の状態への単位時間あたりの散乱確率を示します。したがって、散乱率は

${\displaystyle \lambda (k)=\sum _{k'}S(k,k')}$

これは、第3-3節で説明したように、自由飛行時間と散乱過程を決定するために容易に使用できます。この散乱率は材料のバンド構造に依存することに注意することが重要です（この依存性は行列要素から生じます）。

自由飛行の終わりに、散乱モードと散乱角はランダムに選択されなければならない。散乱メカニズムを決定するためには、すべての散乱率を考慮する必要がある。${\displaystyle \lambda _{1},\lambda _{2},...,\lambda _{n}}$シミュレーションに関連するメカニズム、および散乱時の全散乱率${\displaystyle \lambda _{tot}(t_{sc})=\sum _{i}\lambda _{i}.}$散乱メカニズムを選択すると、0 < r < 1 の一様分布乱数を生成し、以下の規則を参照することになります。

${\displaystyle {\begin{aligned}r&<{\frac {\lambda _{1}}{\lambda _{\mathrm {tot} }}}\rightarrow {\text{scattering-mechanism-}}1\\r&<{\frac {\lambda _{1}+\lambda _{2}}{\lambda _{\mathrm {tot} }}}\rightarrow {\text{scattering-mechanism-}}2\\&{}\ \vdots \\r&<{\frac {\sum _{i=0}^{n}\lambda _{i}}{\lambda _{\mathrm {tot} }}}\rightarrow {\text{scattering-mechanism-}}n\end{aligned}}}$

散乱機構を選択するための計算効率の良いアプローチは、「ボイド」散乱機構を追加することです。${\displaystyle \lambda _{\mathrm {tot} }}$は時間とともに一定に保たれる。粒子がこのメカニズムに従って散乱された場合、散乱後も弾道軌道を維持する。新しい軌道を選択するには、まず散乱後の粒子のエネルギー（または運動量）を導出する必要がある。

${\displaystyle E(k')=E(k)\pm \hbar \omega _{q}\pm \Delta E_{C}\,}$

用語${\displaystyle \hbar \omega _{q}}$フォノン放出または吸収と用語を説明する${\displaystyle \Delta E_{C}}$は、谷間散乱の場合、ゼロではありません。最終エネルギー（およびバンド構造）は、新しい運動量 k' の絶対値を直接与えます。この時点で、散乱粒子の新しい方向（または角度）を選択するだけで済みます。フォノン散乱や放物線状の分散関係などのいくつかの単純なケースでは、散乱角はランダムで、半径 k' の球面上に均等に分布しています。球面座標を使用すると、角度を選択するプロセスは、2 つの角度をランダムに選択することと同等です。${\displaystyle \theta }$そして${\displaystyle \psi }$角度が分布に従う場合${\displaystyle p(\theta ,\psi )}$角度が一様分布の場合、球面上の点を選択する確率は次のようになります。

${\displaystyle p(\theta ,\psi )\,d\theta d\psi ={\frac {\sin \theta \,d\theta \,d\psi }{4\pi }}}$

この場合、2つの変数を分離することが可能です。積分すると${\displaystyle \psi }$そして${\displaystyle \theta }$すると、

${\displaystyle p(\theta )={\frac {\sin \theta }{2}}}$

${\displaystyle p(\psi )={\frac {1}{2\pi }}}$

均一な場合、2 つの球面角は、0 < r ₁、r ₂ < 1 の2 つの乱数を生成することによって選択できます。

${\displaystyle r_{1}=\int _{0}^{\psi }p(\psi ')\,d\psi '={\frac {\psi }{2\pi }}}$

${\displaystyle r_{2}=\int _{0}^{\theta }p(\theta ')\,d\theta '={\frac {1-\cos \theta }{2}}}$

ウィグナー輸送方程式は、ウィグナーに基づく量子補正の基礎となる。

${\displaystyle {\frac {\partial f}{\partial t}}+r\cdot \nabla _{r}f-{\frac {1}{\hbar }}\nabla _{r}V\cdot \nabla _{k}f+\sum _{\alpha =1}^{\infty }{\frac {(-1)^{\alpha +1}}{\hbar 4^{\alpha }(2\alpha +1)!}}\times (\nabla _{r}\nabla _{k})^{2\alpha +1}Vf=\left({\frac {\partial f}{\partial t}}\right)_{c}}$

ここで、kは結晶運動量、V は古典ポテンシャル、右辺の項は衝突の影響、左辺の第 4 項は非局所的な量子力学的効果を表します。左辺の非局所項が空間変化が遅い極限で消えると、標準的なボルツマン輸送方程式が得られます。簡略化された (${\displaystyle \alpha =0}$量子補正されたBTEは次のようになる。

${\displaystyle {\frac {\partial f}{\partial t}}+r\cdot \nabla _{r}f-{\frac {1}{\hbar }}\nabla _{r}V\cdot \nabla _{k}f=\left({\frac {\partial f}{\partial t}}\right)_{c}}$

ここで量子ポテンシャルは項に含まれる${\displaystyle V_{\omega }}$(エラーであるはずです:${\displaystyle V_{\omega }}$（一度も言及されなかった）。

この量子補正法は、1965年にファインマンとヒブスによって開発されました。この方法では、粒子の古典的な経路に沿った量子ゆらぎの経路積分への寄与を計算することによって、有効ポテンシャルが導出されます。この計算は、一次までの試行ポテンシャルを用いた変分法によって行われます。各経路上の平均点における有効古典ポテンシャルは次のようになります。

${\displaystyle V_{\mathrm {eff} }(x)={\frac {1}{{\sqrt {2\pi }}a}}\int _{-\infty }^{\infty }V(x')e^{-{\frac {(x'-x)^{2}}{2a^{2}}}}dx'}$

${\displaystyle a^{2}={\frac {\hbar ^{2}}{12m^{*}k_{B}T}}}$

この手法では、自己無撞着な静電ポテンシャルを入力として、シミュレーションにおいてシュレーディンガー方程式を周期的に解きます。静電ポテンシャルの解に関連する正確なエネルギー準位と波動関数を用いて量子ポテンシャルを計算します。この方法に基づいて得られる量子補正は、次の式で視覚化できます。

${\displaystyle V_{\mathrm {schr} }(z)=-k_{B}T\cdot \log(n_{q}(z))-V_{p}(z)+V_{0}}$

ここで、V _schrは量子補正ポテンシャル、zは界面に垂直な方向、n _qは収束したモンテカルロ濃度に相当するシュレーディンガー方程式からの量子密度、V _pはポアソン解からのポテンシャル、V ₀は半古典的挙動の領域で補正がゼロになるような量子領域から遠く離れた任意の参照ポテンシャルです。上記の量子補正ポテンシャルは計算方法と基本的な仮定が異なりますが、モンテカルロシミュレーションへの組み込みに関してはすべて同じ方法で組み込まれます。