軸対称円錐量子ドット
はじめに
このノートブックではWolfram言語で実装されているように,有限要素法を応用して軸対称形式でシュレーディンガー方程式を解くことを目的とする.ここでの計算は,半導体量子ドット(SQD)に閉じ込められた一粒子状態の解析に焦点を当てたMelnik et al. [Melnik, 2004]の研究に基づいている.ここに提示する方法に従うことで,自身の問題に対して使うことのできる包括的な理解が得られることであろう.
半導体量子ドット(SQD)は2種類以上の半導体材料の結合によって生成された三次元ナノアイランドである.ナノアイランドとそれを取り巻く半導体材料の間にバンド不連続性があるため,半導体量子ドットには伝導電子やホール等の電荷キャリアをド・ブロイ波長に匹敵する大きさの領域に閉じ込めるという特性がある.これが離散エネルギースペクトルになるため,半導体量子ドットは人工原子と呼ばれることもある.このような特徴があるため,SQDはレーザー,太陽電池,単一光子発生器等[Seravalli, 2023] [Sogabe, 2023],さまざまな記述的応用分野で実装されている.
SQDを形成する方法の一つに,Stranski–Krastanov (SK)成長と呼ばれるモードがある.これは分子線エピタキシー法(MBE) [Yu, 2010]という実験的手法に基づいている.このモードでは,等の半導体材料の薄い層を等の半導体基板の上で成長させる.薄い膜が成長し,材料が基板の結晶構造に従うにつれて,ひずみが蓄積し始め,臨界膜厚(通常数ナノメートル)に達すると最終的に自己組織化したアイランドの形成につながる(図1の B) および C) を参照)
図1:SK成長法の説明.図 A) はMBE法で実行される層ごとの成長を示している.図 B) はわずかに変形したによって最初の数層に蓄積しているひずみを示している.図 C) は材料のひずみ緩和による自発的形成を示している.
その後,これらは別の半導体材料(通常基板と同じもの)によって引き続き上を覆うことができ,所望のSQDを得る.特にSK法で成長したSQDは,SQDとほとんど同じ材料からなる「ぬれ層」と呼ばれる薄い膜と結合する.これにもかかわらず,そしてぬれ層の相対的サイズが小さいならば,多くの理論的研究はSQDに閉じ込められた粒子のエネルギースペクトルに対するぬれ層の効果を無視している.それでも閉じ込められた固有状態に対するぬれ層の効果を考慮しないことが適しているかどうかは考慮の余地がある.
さらに成長したSQDの特定の形態は,温度や材料の堆積速度[Granados, 2003]等の形成過程における実験的条件に大きく依存する.例えば,SQDの中には円錐に非常に近い形になるものもあるため,文献ではそのような量子ドットを円錐量子ドット(CQD)と言う.図2を見てみよう.ここでオレンジ色の領域はで更正されるCQDとぬれ層の結合を表し,残りの部分はGaAsから成る.
図2:CQD,ぬれ層,その周りの半導体材料を2つの異なる視点から見た図.ここで ベクトルはシュレーディンガー方程式の電子の位置を表す.
先に述べた通り,このドキュメントの中の計算はMelnik et al. [Melnik, 2004]の研究に基づいている.これはとの結合から成るCQD/ぬれ層構造に閉じ込められた一粒子状態を研究したものである.ここでは特に,時間非依存のシュレーディンガー方程式およびそれに対応する固有値問題の解き方を示す.その結果として,幾何学的特性に閉じ込められた一粒子状態に対するぬれ層の影響も調べることができる.Wolfram言語に実装されている有限要素法を使う利点を示すことができるであろう.
理論
有効質量と包絡関数近似(詳細はこちら)を考えると,電子の波動関数についての時間非依存のシュレーディンガー方程式を方程式1として書くことができる.
ここで は閉込めポテンシャルである.これはCQDとぬれ層の領域内では値を取り,それ以外では有限の一定値を取る.さらに有効質量関数 はそれぞれの材料に関連付けられた個別の有効質量を与える区分関数である.
方程式2の運動項(拡散項)の書き方に注目する.ハミルトニアンをこのように書くことを通常BenDaniel-Dukeハミルトニアンと言う.これは各領域における有効質量間の不一致に対してより正確なモデルを与える[BenDaniel-Duke, 1966].
図3:Melnikの研究によるCQD領域の断面図[Melnik, 2004].数字の1から6はそれぞれの境界面を表す.
図3では,問題の全体の領域が平面の横断面として描かれている.ここで数字のからは異なる境界面を表す.面とで囲まれた体積は,から成るCQDとぬれ層の結合に該当する.領域の残りの部分はでできている.
特に興味があるのは束縛状態,つまりCQDとぬれ層から十分遠いところで減衰する状態だけである.このため,面とで確実にとなるディリクレ境界条件条件を使う.対照的に,ぬれ層全体で分散しCQDの閉込めを超えて拡張する状態もあるので,面とにはディリクレ境界条件を使わないようにする.その代りにノイマン境界条件()を使う.その上,境界とでBenDaniel-Duke境界条件を強制する必要がある.つまり関係 がこれらの内部境界で真とならなければならないのである[BenDaniel-Duke, 1966].
一方,形状の円筒対称性を利用し,閉込めポテンシャルと有効質量関数の両方が方位角 とは関係ないと仮定すると,問題の次元をDからDに減らすことができる.
関数が に依存しないと仮定すると,円筒座標の微分演算子の定義を適用することで以下が得られる:
CQDの 軸についての円筒対称性があるとすると,仮説 が提案できる.以下が導ける:
ここで方程式全体を で割ると,問題を2つの方程式に分けることができる.並べ替えて簡約すると以下が得られる:
ここで軌道角運動量または角運動量量子数 は周期性のため整数である.
このアプローチにより,図4で示すように方程式5についてD幾何で解き,エネルギー固有値 とD波動関数 を得ることができる.したがって,全は導関数は方程式3で得られる.
形状
この系をモデル化する最初のステップはパラメータの定義である.長さはすべてナノメートルとする.
ここでWLWidthはぬれ層の厚さである.dotRadiusとdotHeightはそれぞれCQDの底面半径と高さである.domainHeightとdomainRadiusはそれぞれこれから解く方程式5のD領域全体の高さと幅である.
ここで領域を定義する.特にCQDとぬれ層が結合した領域と残りの領域を区別したいので,別々の領域を定義する.
CQDに対応する領域ではTriangle関数を使い,ぬれ層の領域ではRectangle関数を使う.それからRegionUnionを使って一つのCQD/ぬれ層領域を作成する.の値がぬれ層の幅の中点と整列することに注目する.
解の品質をよりよく制御するために,まず有限要素法に使うメッシュを生成する.これによってメッシュをより簡単に可視化し,必要に応じて調整することができる.
RegionBoundary関数を使ってCQD/ぬれ層領域の境界を抽出し,それを組み合せて一つの境界メッシュにする.
関心があるのは束縛状態なので,波動関数が大きい値を取りそうなところ(この場合CQD/ぬれ層領域内部)でメッシュを調整したい.そのためにToElementMesh関数の"RegionMarker"オプションを使う.こうすることで,要素マーカーの領域に対する調整パラメータを定義することができる.
ハミルトニアン
次のステップではPDE演算子を作成する必要がある.このためにポテンシャルと有効質量関数を定義する.
閉じ込めポテンシャルは,それぞれの材料のエネルギーギャップ間の違いによって引き起こされる.InAsとGaAsの結合の場合,ポテンシャル障壁は値 V0=によって近似することができる[Melnik, 2004].
閉じ込めポテンシャルを定義するために,先のメッシュの定義で使った要素マーカーを利用する.
InAsとGaAsの有効質量はそれぞれ ,である[Yu, 2010].ここで は通常の電子質量である.
ここでは方程式6を解くことが目標である.したがって対応するPDE演算子を生成するためにSchrodingerPDEComponent関数を使う.
パラメータ"RegionSymmetry" -> "Axisymmetric"と"AzimuthalQuantumNumber" -> lを使っていることを特筆すべきであろう.この2つのパラメータによって,SchrodingerPDEComponentは方程式7の項 および に対応する式を生成する.またすべての形状パラメータはナノメートルで定義されているので,パラメータ"ScaleUnits" -> {"Meters" -> "Nanometers"}を使う.つまりモデルパラメータが変換され,長さの単位の"Meters"から"Nanometers"に変わったということである.このことでPDE演算子がSchrodingerPDEComponent関数から正しく生成されていることが分かる.
境界条件と固有問題の解
先に説明したように,CQD/WL領域の束縛条件だけが見付かるようにディリクレ条件を設定する.
またノイマン境界条件を図2の境界に適用しなければならない.この場合はデフォルトで適用されているので,この境界条件を明示的に定義する必要はない(NeumannValueを参照のこと).また数量 はデフォルトでは内部境界の間で連続である.つまりBenDaniel-Duke境界条件も正しく適用されているということである.
ここでNDEigensystemを使って固有値問題を解くことができる.のときの最初のつの低い固有状態から始める.
まずエネルギーを正しいエネルギー単位に変換する必要がある.モデルのパラメータの長さの単位がとすると,エネルギー単位はがに変わることを除いて"SIBase"単位である.
さらに変数funs[[i]]に対応した方程式8から9までの波動関数に関連付けられた確率密度を可視化することができる.
プロットからいくつかの点についてコメントすることができる.まず,基底状態では確率密度は予想通り主にCQDの中心に閉じ込められている.一方,他の状態では確率密度は基底状態よりもぬれ層で拡散している.
軌道角運動量の影響
軌道角運動量の影響を調べるために,まず の値をと設定し,前と同じ方法で固有値問題を解く. にその値を選んだ理由は後で分かる.
の場合,最後の固有値はポテンシャル障壁(V0=)より大きいので,非束縛状態である.これは確率密度をプロットすることで確認することができる.
一つ気を付けなければならないのは, を非零の値に設定することによって,波動関数は ではなく になるということである.つまり,i 番目の状態が式funs[[i]]Exp[lϕ]で表されるので波動関数は複素数値になるのである.それにも関わらず確率密度はまだ なので i 番目の確率密度はfuns[[i]]2で表される.
ここで軌道角運動量のおもしろい一面に気が付く.まず番目の固有状態は予想通り非束縛である.さらにどの固有状態の場合でも における確率密度はである.特に方程式10の項 は反発項であり,この動作を説明している. がに近付くと,この項は無限に近付くため,波動関数および確率密度が においてに減衰するよう強制する.また,それぞれの確率密度で分かるように,この影響は l=0では存在しない.この動作は の非零の値で見られるので,の値を選ぶ必要があったのである.
一方,DensityPlot3Dを使って完全な波動関数を可視化することができる.そのためには, に非零の値を導入したため波動関数が複素数値になっていることを考慮に入れなければならない.それゆえ,虚部と実部を別々にプロットする必要がある.
また, を変更することでエネルギー固有値にどのような影響があるかを調べることもできる.そのために,ヘルパー関数を定義してそれぞれの の値に対する固有値問題を解く.
をからの範囲で変化させたときのエネルギー固有値をプロットすると, の絶対値が大きければ大きいほどエネルギー固有値も大きいことが分かる.これは が角運動量量子数と考えると説明できる.CQD/ぬれ層構造の対称軸である 軸を中心とする角運動量に関連付けられた量子機械演算子は方程式11で書くことができる.
方程式12の波動関数が方程式13の演算子 の固有関数であることは簡単に分かる.したがって,量子力学によると,対応する固有値 は 軸についての角運動量の可能な値を表す.その結果, の大きい値は 軸に対して大きい角運動量を持つ状態に直接対応する.同様にこれは の値を増やしたときにエネルギー固有値が大きくなることについての物理的説明になっている.
一方,赤い横線が障壁の高さ V0=を表す上のプロットから分かるように, の値が大きいほどCQD/WL領域に束縛される状態は少なくなる.これはよりエネルギーのある状態なので,非束縛状態が多くなることが予測される.
ぬれ層のサイズの影響
ぬれ層の大きさの変化の影響を調べるために,ぬれ層の幅以外の形状パラメータを固定しておく.
簡単さのため,そして を変化されることの影響は上ですでに調べたため,における軸対称量子数に対するぬれ層の幅の関数として固有値を計算する. とぬれ層の幅の両方を変化させたときの影響を調べたい場合は,固有値をWLWidth パラメータの関数として計算するように以下のヘルパー関数の定義を変更するとよい.
さまざまな形状設定でのエネルギー固有値を計算するためには,新しいメッシュを作成し,閉じ込めポテンシャルと有効質量関数を再定義する必要がある.
最後のプロットでは,ぬれ層の幅が小さくなると,エネルギー固有値の縮退が見られる.例えば,パラメータWLWidthがからまでの値では状態およびが縮退している.この縮退はぬれ層の幅がもっと大きくなるとなくなる.さらにぬれ層の幅の値がおよびでは,WLWidthの低い値の場合と異なり,状態が束縛状態になる.
概して,ぬれ層の幅が増加するにつれて固有値は単調に減少する.電子が閉じ込められる全体積が大きくなるので,想定できることである.ハイゼンベルクの不確定性原理を踏まえると,ぬれ層の幅の値が小さくなると波動関数がより充填されるためその位置の不確定性は最低になる.つまり,運動量における不確実性は最高になるということである.換言すると,粒子がより多く閉じ込められるとその平均平方運動量が大きくなり運動エネルギーも大きくなる.このような理由で,ぬれ層の幅の値が低いほどエネルギー固有値が大きくなるのである.また逆の影響もある.ぬれ層が大きくなるにつれ波動関数はより拡散し,その位置の不確実性は増加する.これが運動量の不確実性を減少させる.平均平方運動量は低くなり,これは粒子に対する運動エネルギーが低いことを意味する.
これを実証するために,ぬれ層がのときの結果とのときの結果の束縛状態の確率密度を比較することができる.ぬれ層の幅がのときの束縛状態は最初のつの状態だけなので,それだけを選ぶ.
さまざまな形状パラメータに相当する波動関数をプロットするために,各解からメッシュが抽出されていることに注目のこと.
上では最初のつの束縛状態の確率密度がプロットされている.左の列と右の列はそれぞれWLWidthパラメータの値がとに対応している.ぬれ層の幅の値が小さいほど確率密度が圧縮され局所化が大きくなることが分かる.ぬれ層が大きい場合の波動関数と比べて,低位状態では特にそうである.
まとめ
ここまでぬれ層と結合したCQD内に閉じ込められた一電子系の固有状態対する軌道角運動量とぬれ層の厚さの影響を調べてきた.軌道角運動量の導入により波動関数が対称軸で消失することが分かった.さらに軌道角運動量の値が増加するとエネルギー固有値も増加する.
またぬれ層の幅が大きくなるとエネルギー固有値は減少した.これは先に説明した,想定通りの動作である.この動作は,ぬれ層の大きさがエネルギー固有値に大きな影響を及ぼしており,モデルにぬれ層を組み込むことの重要性を強調している.
ここではWolfram言語を使って同様の問題を解こうとしている方のために主なメソッドを説明した.
参考文献
1. Seravalli, L. (2023). Metamorphic InAs/InGaAs quantum dots for optoelectronic devices: A review. Microelectronic Engineering, 276, 111996. https://doi.org/10.1016/j.mee.2023.111996
2. Sogabe, T., Shoji, Y., Miyashita, N., Farrell, D. J., Shiba, K., Hong, H.-F., & Okada, Y. (2023). High-efficiency InAs/GaAs quantum dot intermediate band solar cell achieved through current constraint engineering. Next Materials, 1(2), 100013. https://doi.org/10.1016/j.nxmate.2023.100013
3. Yu, P. Y., & Cardona, M. (2010). Fundamentals of Semiconductors (pp. 8–13). https://doi.org/10.1007/978-3-642-00710-1
4. Granados, D., Garcia, J. M. (2003). Customized nanostructures MBE growth: From quantum dots to quantum rings. Journal of Crystal Growth, 251(1–4), 213–217. https://doi.org/10.1016/S0022-0248(02)02512-5
5. Melnik, R. V. N., Willatzen, M. (2004). Bandstructures of conical quantum dots with wetting layers. Nanotechnology, 15 (1), 1–8. https://doi.org/10.1088/0957-4484/15/1/001
6. BenDaniel, D. J., Duke, C. B. (1966). Space-Charge Effects on Electron Tunneling. Physical Review, 152 (2), 683–692. https://doi.org/10.1103/PhysRev.152.683