シリコン／SiGe量子井戸：六バンドルッティンガー・コーン模型によるホール輸送シミュレーションの解説
1. はじめに
当社のシリコンΔバレー・チュートリアルでは、Si/SiGe量子井戸内の電子挙動について解説した。本チュートリアルでは、その電子構造のもう一方の側面、すなわち価電子帯内のホール挙動を取り扱う。層構造自体は類似しているものの、基礎となる物理現象は大きく異なる。個々のシリコン伝導バレー付近の電子は異方性を持つ有効質量を用いて記述可能だが、価電子帯の頂点付近には互いに相互作用する複数の状態が存在する。したがって、単一のホール分散パラボラモデルでは重要な情報が見落とされることになる。
ここではOghmaNanoの6バンドLuttinger–Kohnモデル（LK6）を用いて、厚さ10nmのシリコン量子井戸内におけるホールのサブバンド分布を計算する。さらにバンドエネルギーの算出にとどまらず、音響フォノン散乱やホールの占有状態、そしてそれらに基づく面内移動度も評価する。構造を一から作成することなく、提供されているデフォルトの例をそのまま利用することも可能である。
2. 六バンドモデルの例を開く
「ファイル」→「新規シミュレーション」を選択します。サンプルブラウザ内で「GaAsデモ」(??)をダブルクリックします。このフォルダにはシリコン関連のサンプルを含む半導体量子井戸用テンプレートも含まれています。次に表示されるウィンドウで「ひずみを伴うSi/SiGe価電子帯量子井戸（LK6）」(??)をダブルクリックします。必ず「価電子帯」用のテンプレートを選択してください。隣接する「Delta-Valley」テンプレートはシリコン中の伝導電子を計算するためのものです。
3. Si/SiGe構造および材料パラメータの確認
メインウィンドウにはサンプルの三次元表示が示されている（??）。電荷キャリアは、我々がz方向と呼ぶ層方向に閉じ込められる一方で、x-y平面内では自由に移動できる。［レイヤーエディター］をクリックすると、各層の実際の厚みや中央層が量子井戸として設定されている様子を確認できる。
| 層 | 厚み | 役割 |
|---|---|---|
| シリコンゲルマニウム | 20ナノメートル | 周囲層 |
| シリコン（量子井戸） | 10ナノメートル | 量子閉じ込めされた価電子状態 |
| シリコンゲルマニウム | 20ナノメートル | 周囲層 |
??に示されている光学材料の設定欄は、電子材料の定義内容を確認するための代替手段ではありません。バンド端のオフセットや合金組成、ひずみ設定などが実際の量子閉じ込めポテンシャルに影響を与えるため、初回のシミュレーションでは例文で示されている設定内容をそのまま使用してください。
電気的パラメータ編集画面を開き、Si (QW)を選択してください（??）。スクリーンショットにはLuttingerパラメータγ1およびその文献情報が示されている。単一バンドの電子モデルとは異なり、価電子帯の計算では重ホール、軽ホール、スプリットオフ帯を記述するために3つのLuttingerパラメータが用いられる。さらにスピン軌道相互作用による分裂の考慮が必要であり、ひずみが有効な場合には適切な変形ポテンシャルパラメータも必要となる。
4. なぜ価電子帯の記述には6つの成分が必要なのか？
シリコンの価電子帯頂点付近には、重ホール（HH）および軽ホール（LH）という2つのホール状態群が存在する。さらにスピン軌道相互作用により、よりエネルギーの低いスプリットオフ帯（SO）も形成される。各状態群にはスピンに関連する2つの基底状態が存在するため、LK6モデルでは合計6つの成分が得られる。これらは互いに結合した波動関数の基底成分であり、必ずしも独立して運動する6つの放物線型バンドを意味するわけではない。
量子井戸の状態は、6成分からなる列ベクトルとして模式的に表現できる。
\[ \Psi(z)=\begin{pmatrix} F_{\mathrm{HH}+}(z) & F_{\mathrm{LH}+}(z) & F_{\mathrm{SO}+}(z) & F_{\mathrm{HH}-}(z) & F_{\mathrm{LH}-}(z) & F_{\mathrm{SO}-}(z) \end{pmatrix}^{\!T}. \]
ここで用いられているラベルは、基礎となる角運動量基底の主要な特性を示している。面内運動量が変化するにつれて各成分は一般に混合されるため、初期状態が重ホール的と記述された状態であっても、より大きな波数ベクトルにおいては必ずしも純粋な重ホール状態のままであるとは限らない。コンピュータによって解かれる問題は以下の形式となる。
\[ \left[H_{\mathrm{LK}}(k_x,k_y,-i\partial_z) +H_{\mathrm{BP}}(\varepsilon)+V_v(z)I_6\right]\Psi =E(k_x,k_y)\Psi. \]
Luttinger–Kohn項は価電子帯状態の運動エネルギーおよび状態間の混合を記述し、Bir–Pikus項はそれらがひずみに対してどのように応答するかを示す。また\(V_v(z)\)は価電子帯の量子閉じ込めポテンシャル形状を表す。本モデルでは、重ホール、軽ホール、スプリットオフ帯それぞれに独立した質量定数を割り当てるのではなく、基底ベクトル成分間の結合を維持している。本チュートリアルにおいてはこの行列を手動で構築する必要はない。詳細なハミルトニアンについてはLK6理論ページを参照されたい。
比較のため、単純化されたバルク近似を用いても、なぜ3つのルッティンガーパラメータが重要であるかが明らかになる。すなわち、立方晶の主軸方向において、重ホールおよび軽ホールの分散曲線の曲率は、これらのパラメータの異なる組み合わせと関連しているのである。
\[ \frac{m_0}{m_{\mathrm{HH}}^{*}}\simeq\gamma_1-2\gamma_2, \qquad \frac{m_0}{m_{\mathrm{LH}}^{*}}\simeq\gamma_1+2\gamma_2. \]
これらは説明的な極限関係として有用なものに過ぎず、完全な六バンド量子井戸計算を二つの別々の放物線で置き換えるよう指示するものではない。閉じ込め構造においては、井戸幅、スピン軌道相互作用、ひずみ、および面内運動量が実際の固有値ならびにその状態特性の変化を決定するのである。
5. 計算の実行
ファイルリボンに戻り、「シミュレーション実行」をクリックしてください。OghmaNanoはまず、一連の面内波数ベクトルに対して、量子閉じ込めされたLK6価電子帯状態を求めます。計算の進行状況はターミナルタブ（??）で確認可能です。表示されるエネルギー値（target）は内部固有値計算の目標値であり、出力される収束メッセージは要求された固有状態が数値的に求められたことを示しています。
ソルバーの計算が完了するのを待ってください。散乱過程の計算にはより多くの時間が必要となります。なぜなら、単にサブバンドのエネルギー値を求めるだけでなく、占有されている初期電子状態から利用可能な最終電子状態への遷移も評価しなければならないからです。
6. 計算された価電子帯の分散関係の読み取り方
「出力」タブを開いてください。生成されるファイルには、バンド関連のプロットのほか、散乱および移動度に関するフォルダーが含まれています（??）。価電子帯の結果ファイルを開くと、??が得られます。
\(k_{\parallel}=0\)においては、量子閉じ込め効果によって各サブバンドの開始エネルギーが離散化される。有限な\(k_{\parallel}\)領域に移ると、示されている電子エネルギーの定義に従い各エネルギー曲線は下向きに湾曲し、その曲率や相互間隔もやや異なっている。近接して存在する線群は、時間反転対称性を有する計算におけるスピン関連の（クラマース）縮退ペアに対応しており、これらを自動的に異なる種類のホールとして扱うべきではない。
このプロットを理解するための有効な方法は、ホールとは価電子帯における電子の欠落状態であることを念頭に置くことです。したがって、ホールのエネルギーは選定された価電子帯の基準値を基準として表現することができます。
\[ E_h(\mathbf{k}_{\parallel})=E_{v,\mathrm{ref}}-E_v(\mathbf{k}_{\parallel}). \]
したがって、最大値から下向きに湾曲する価電子のエネルギー曲線は、ホールの励起エネルギーが増加していることを示す。各曲線が持つ重ホール・軽ホール・スプリットオフ帯の特性は、単にその色や曲率からではなく、6成分からなる固有ベクトルによって決定される。
この結果をシリコンΔバレー・チュートリアルで説明されている伝導帯の分散と比較してみる価値がある。そちらでは異なるバレーの向きを異方性を持つ電子質量を用いて扱っているが、本モデルでは多バンド価電子ハミルトニアンにおいてホール状態同士が直接結合しているのである。
7. 音響フォノン散乱とホールの寿命
バンド構造からはどのようなホール状態が存在しうるかは分かりますが、それだけではホールがどの程度速く運動量を失うかを予測することはできません。本例においてOghmaNanoは、計算によって得られた閉じ込め状態を用いて音響フォノン散乱の解析を行っています。音響フォノンは結晶内に微小なひずみ場を生じさせ、価電子帯を摂動させることでホールを初期状態\(i\)から最終状態\(f\)へと移動させる可能性があります。
si_lk6_acoustic_rate_*.csvやsi_lk6_acoustic_tau_*.csvをはじめ、輸送特性に関するデータセットも含まれている。「出力」タブ内の「散乱」フォルダーを開いてください（??）。番号付きのファイル群が計算結果の状態を表しています。例えば、si_lk6_acoustic_tau_02.csvのスクリーンショットでは、音響フォノン散乱によるホールの寿命と面内波数ベクトルとの関係がプロットされています。
物理的な出発点としてはフェルミの黄金律が挙げられる。
\[ W_{i\rightarrow f}=\frac{2\pi}{\hbar} \left|\langle\Psi_f|H_{\mathrm{e-ph}}|\Psi_i\rangle\right|^2 \delta(E_f-E_i\pm\hbar\omega_{\mathbf q}). \]
この行列要素は、格子振動が2つのホール状態をどの程度強く結合させるかを示しており、デルタ関数はフォノンの吸収または放出におけるエネルギー保存則を表している。6バンド価電子帯モデルにおいては、この結合の強さは複数成分からなるホール波動関数に依存するため、許容される遷移は重ホール・軽ホール・スプリットオフ帯間の混合ならびに利用可能な終状態のエネルギー値にも影響される。ここで出力される結果は、本例で用いられている音響フォノン近似に関するものである。
si_lk6_acoustic_tau_02.csvに示す音響フォノン散乱によるホールの寿命値である。この値が波数ベクトルによって変化するのは、利用可能な終状態およびそれらとの結合の度合いに起因する。許容される遷移をすべて合計することで、全事象発生率 \(\Gamma_i\) が得られ、それに対応する寿命も求められる。
\[ \Gamma_i=\sum_f W_{i\rightarrow f},\qquad \tau_i=\frac{1}{\Gamma_i}\quad(\Gamma_i\gt 0). \]
??に示されているホールの寿命値は一定ではない。波数ベクトルが結合した価電子帯分散曲線上の各領域を移動するにつれてこの値は変化する。最初のサンプリング点において得られるゼロ値は、物理的には\(k_{\parallel}=0\)において散乱寿命が完全にゼロであることを示すものではなく、境界条件や出力形式上の表記、あるいは検証が必要な値として扱うべきである。
輸送現象においてはさらに重要な区別があり、それはすべての衝突が電流を同等に緩和させるわけではないという点である。小角度散乱は微視的な状態を変化させる可能性があるが、運動の全体的な方向性には比較的影響を与えない。このため、輸送重み付けされた散乱率はしばしば次のように記述される。
\[ \Gamma_{\mathrm{tr},i}=\sum_f W_{i\rightarrow f} \left(1-\cos\theta_{if}\right),\qquad \tau_{\mathrm{tr},i}=\frac{1}{\Gamma_{\mathrm{tr},i}}. \]
この角度に関する表現は運動量緩和の原理を示しており、異方性を持つ混合バンドの場合には輸送計算全体において適切な速度値および重み付けを用いる必要がある。このため、出力結果では音響フォノン散乱の確率や寿命が輸送に関連する物理量から分離されて示されており、単一の数値があらゆる目的に用いられることはないのである。
8. ホール密度およびフェルミ準位から移動度を求める方法
測定可能な輸送特性を得るためには、散乱に関する情報と、各サブバンドの各部分にどれだけのホールが存在するかを組み合わせる必要がある。そこで本チュートリアルでは、単位面積当たりのホール密度を示す\(p_{2D}\)の値を変化させており、その単位はm⁻²である。
si_lk6_hole_fermi_level.csv におけるホールのフェルミ準位。
si_lk6_hole_mobility.csvから算出されたホールの移動度。フェルミ準位は、どの電子状態が占有されているかを決定する。価電子のエネルギーが\(E_v\)である場合、対応するホール状態の平衡時の占有数は次式で表される
\[ f_h(E_v)=1-f_e(E_v) =\frac{1}{1+\exp\!\left[(E_F-E_v)/(k_B T)\right]}. \]
シートホール密度を高めると、化学ポテンシャルが価電子帯状態に対して移動することになり、これが左側の図で示されている。具体的なエネルギー値はモデルで用いられる基準値に依存するが、エネルギー掃引の方向およびそれに伴うホールの分布状態が最も重要な特徴となる。
OghmaNanoは、占有状態におけるキャリアの速度および輸送緩和時間に重み付けを行い、面内のシート伝導率を求める。 \(\sigma_{2D}\)これに対応するホールの移動度は
\[ \mu_h=\frac{\sigma_{2D}}{q\,p_{2D}}, \qquad \sigma_{2D}\sim q^2\sum_n\int \frac{d^2 k_{\parallel}}{(2\pi)^2} \,v_{n,x}^{2}\,\tau_{\mathrm{tr},n} \left(-\frac{\partial f_e}{\partial E}\right). \]
ここで\(v_{n,x}\)はサブバンド\(n\)における面内群速度を表し、積分記号は異なる波数ベクトルを持つ状態が電流生成に異なる寄与をすることを示している。移動度のグラフ(??)では、示された密度範囲内（およそ0.165から0.151 m2 V−1 s−1）で移動度がわずかに減少している。これは設定済みの例およびそこに含まれる散乱物理モデルに基づく結果であり、普遍的なシリコンのホール移動度則やあらゆる実験的損失機構を考慮した予測値ではない。
9. 自分で試してみよう
💡 推奨される実験例：デフォルトのサンプルを保存した上で、一度に1つのパラメータのみを変更してみてください。
- Si量子井戸の幅を10 nmから変更し、計算を再実行してください。その後、閉じ込められた価電子帯のエネルギー準位を比較してみましょう。一般的に井戸幅が狭くなるほど量子閉じ込めの効果は強くなりますが、重ホール・軽ホール・スプリットオフ帯の各状態が相互に影響し合うため、その応答は単一質量モデルよりも複雑になります。
- ホールモデルと電子モデルを比較するには、同様の層構造を用いて付属のシリコン Δバレー例を実行してください。価電子帯用のソルバーでは6成分からなるホール状態が連成して記述されるのに対し、伝導帯の計算ではシリコンのバレー方向が個別に区別される点に注目してください。
si_lk6_acoustic_tau_00.csvといった別の音響フォノン散乱寿命ファイルを確認し、si_lk6_acoustic_tau_02.csvと比較してみましょう。異なるサブバンドにおいては、利用可能な散乱チャネルが同一である必要はありません。- 密度スキャンの結果を確認してください。単一の有効寿命値で両方の曲線が予測できると仮定するのではなく、フェルミ準位と移動度の変化を比較する必要があります。
これにて、六バンドモデルを用いた量子井戸に関する入門シリーズは完了となる。完全な数学的定式化および各物質パラメータの意味については、六バンドルッティンガー・コーン理論の章およびバー・ピクス歪み理論の章に記述されている。