Розуміння теорії k·p: від параболічних зон до багатозонних моделей.

Навіщо нам потрібна ще одна модель зонної структури?

У нашому вступі до зон параболічної квантової ями ми представили енергію електрона у вигляді простої параболи. Її кривизна визначалася ефективною масою, що дозволило нам перетворити кожен локалізований квантово-ямовий стан на двовимірну електронну підзону. Це дуже корисне наближення, особливо в області, близькій до нижньої межі простої зони провідності.

Однак, реальний напівпровідник містить декілька електронних зон. Ці зони можуть розташовуватися близько одна до одної, відхилятися від параболічної форми та взаємодіяти між собою. Це особливо важливо для валентних зон, де важкі, легкі та відщеплені стани не завжди можна розглядати як незалежні частинки з єдиною масою.

Метод k·p надає можливість описувати ці ефекти без повторного обчислення повної електронної структури на атомному рівні для кожного хвильового вектора. Ми починаємо з електронних станів у одній добре вивченій точці зонної структури кристала та обчислюємо, як вони змінюються під час руху хвильового вектора від цієї точки. Найважливіше те, що ми дозволяємо різним електронним зонам взаємодіяти між собою.

Мета цієї сторінки – пояснити, звідки беруться рівняння k·p, чому ефективна маса виникає природним чином і як система пов’язаних рівнянь перетворюється на багатосмуговий гамільтоніан, який використовується в програмі для розв’язання задач, пов’язаних із квантовими ямами.

Почнемо з електрона в кристалі.

Для опису руху електрона в кристалі напівпровідника ми почнемо з добре відомого рівняння Шредінгера, яке не залежить від часу:

\[ \hat H_0\Psi(\mathbf r)=E\Psi(\mathbf r), \qquad \hat H_0=\frac{\hat{\mathbf p}^{\,2}}{2m_0}+V_{\mathrm{crystal}}(\mathbf r). \]

Тут \(m_0\) – це звичайна (маса вільних електронів), \(\hat{\mathbf p}=-i\hbar\nabla\) – оператор імпульсу, а \(V_{\mathrm{crystal}}\) – потенціал, який відчуває електрон через періодичне розташування атомів. Ми свідомо використали \(m_0\), а не ефективну масу напівпровідника: на цьому етапі вплив кристалічної решітки ще не враховано в наближеній масі.

Важливою властивістю ідеального кристала є те, що його потенціал повторюється від однієї елементарної комірки до наступної. Якщо \(\mathbf R\) описує трансляцію кристалічної ґратки, то:

\[ V_{\mathrm{crystal}}(\mathbf r+\mathbf R) =V_{\mathrm{crystal}}(\mathbf r). \]

Ми могли б спробувати розв’язати рівняння Шредінгера безпосередньо для всієї цієї періодичної потенціальної області, але теорема Блоха пропонує набагато зручніший спосіб представлення її розв’язків.

Теорема Блоха: відокремлення хвилі від кристалічної решітки

Теорема Блоха стверджує, що електронний стан у періодичному кристалі можна представити у вигляді добутку площинної хвилі та функції, яка повторюється відповідно до кристалічної ґратки:

\[ \boxed{ \Psi_{n\mathbf k}(\mathbf r) =e^{i\mathbf k\cdot\mathbf r}\,u_{n\mathbf k}(\mathbf r) } \]

Обидві частини виконують різні функції. Фактор площинної хвилі описує, як загальна фаза змінюється під час проходження через кристал. Функція \(u_{n\mathbf k}(\mathbf r)\) описує деталі хвилі всередині кожної елементарної комірки та має таку ж періодичність, як і кристал:

\[ u_{n\mathbf k}(\mathbf r+\mathbf R) =u_{n\mathbf k}(\mathbf r). \]

Позначка \(n\) вказує на те, яку електронну зону ми розглядаємо, а \(\mathbf k=(k_x,k_y,k_z)\) – це хвильовий вектор кристала. Величину \(\hbar\mathbf k\) називають кристалічним імпульсом. Вона тісно пов’язана з імпульсом, який ми використовували в нашій попередній аналогії з крикетним м’ячем, хоча ці дві величини не слід плутати в періодичному кристалі.

Це перше важливе спрощення: замість того, щоб описувати нову хвильову функцію атомного масштабу для всієї кристалічної решітки для кожного \(\mathbf k\), ми можемо працювати з періодичними функціями, що повторюються в межах елементарної комірки \(u\). Метод k·p досліджує, як ці функції та їхні енергії змінюються, коли ми переміщуємося на невелику відстань у просторі хвильових векторів.

Звідки береться член k·p у формулі?

Підставимо вираз Блоха в рівняння Шредінгера. Коли оператор імпульсу діє на добуток площинної хвилі та періодичної функції, визначеної в комірці, результатом буде:

\[ \hat{\mathbf p}\left(e^{i\mathbf k\cdot\mathbf r}u\right) =e^{i\mathbf k\cdot\mathbf r} \left(\hat{\mathbf p}+\hbar\mathbf k\right)u. \]

Іншими словами, дія на блохівську хвилю еквівалентна заміні оператора імпульсу, що діє на \(u\), на \(\hat{\mathbf p}+\hbar\mathbf k\). Застосувавши це до члена, що описує кінетичну енергію, і скоротивши спільний множник у вигляді площинної хвилі, отримаємо:

\[ \left[ \frac{(\hat{\mathbf p}+\hbar\mathbf k)^2}{2m_0} +V_{\mathrm{crystal}}(\mathbf r) \right]u_{n\mathbf k} =E_n(\mathbf k)u_{n\mathbf k}. \]

Розширення квадрата є вирішальним кроком. Оскільки \(\mathbf k\) є параметром, а не просторовою похідною, ми отримуємо:

\[ \boxed{ \left[ \hat H_0 +\frac{\hbar}{m_0}\mathbf k\cdot\hat{\mathbf p} +\frac{\hbar^2 k^2}{2m_0} \right]u_{n\mathbf k} =E_n(\mathbf k)u_{n\mathbf k} } \]

Перший доданок – це наш початковий кристалічний гамільтоніан. Останній доданок – це проста квадратична складова енергії. Середній доданок, \(\hbar\mathbf k\cdot\hat{\mathbf p}/m_0\), є основою для назви методу: k·p.

З фізичної точки зору, це показує, що відхилення від обраного опорного хвильового вектора може призвести до змішування електронних станів за допомогою оператора імпульсу. Щоб зрозуміти, що це означає, нам потрібно визначити, між якими станами може відбуватися це змішування.

Вибір орієнтиру та основи.

Структура зон провідності напівпровідників часто є найбільш цікавою поблизу мінімуму або максимуму зони. Для напівпровідників із прямою зоною, таких як GaAs, мінімум зони провідності та максимум валентної зони знаходяться в центрі зони, який називається точкою Гамма (\(\Gamma\)). У подальшому поясненні ми, отже, розглядатимемо область навколо точки \(\mathbf k=0\). Ту саму загальну ідею можна застосувати до іншої опорної точки, коли відповідний екстремум зони знаходиться в іншій частині зони Бріллюена.

У точці \(\mathbf k=0\) член k·p зникає. Припустимо, що ми вже знаємо декілька кристалічних станів, які періодично повторюються в межах елементарної комірки:

\[ \hat H_0|u_{n0}\rangle=E_n(0)|u_{n0}\rangle. \]

Ці відомі стани є основою для наших розрахунків. Стан, що відповідає близькому хвильовому вектору, можна описати як їхню комбінацію:

\[ |u_{\mathbf k}\rangle =c_1|u_{10}\rangle+c_2|u_{20}\rangle+\cdots. \]

Коефіцієнти \(c_1,c_2,\ldots\) показують, яка частина кожного початкового стану впливає на новий стан. Саме їх потрібно точно визначити в наших розрахунках. Ми робимо це, підставляючи розкладання в рівняння k·p і проєктуючи його на кожен базовий стан.

У результаті отримуємо матричне рівняння. Його елементи:

\[ \boxed{ H_{nm}(\mathbf k) =\left[E_n(0)+\frac{\hbar^2k^2}{2m_0}\right]\delta_{nm} +\frac{\hbar}{m_0}\mathbf k\cdot\mathbf p_{nm} } \]

\[ \mathbf p_{nm} =\langle u_{n0}|\hat{\mathbf p}|u_{m0}\rangle. \]

Символ \(\delta_{nm}\) дорівнює одиниці, коли дві смуги мають однакові позначення, і нулю в інших випадках. Отже, діагональ цієї матриці містить еталонні значення енергії смуг і квадратичні члени, що описують вільні електрони. Її недіагональні елементи містять інформацію про взаємодію імпульсів між різними смугами.

Ці елементи матриці імпульсу є ключовими для теорії k·p. Якщо значення параметра взаємодії дорівнює нулю, відповідні стани не взаємодіють через цей член. Якщо значення параметра взаємодії є значним, їхні енергії та хвильові функції можуть суттєво змінюватися при віддаленні від опорної точки.

Практичний приклад: що відбувається, коли дві музичні групи взаємодіють?

Перш ніж розглядати вісім або десять зон, корисно подивитися, як працює модель із двома зонами. Уявімо, що ми зберігаємо лише один стан провідної зони та один стан валентної зони, з еталонними енергіями \(E_c\) і \(E_v\). Їхня різниця в еталонній точці є шириною забороненої зони, \(E_g=E_c-E_v\).

Для наочності розглянемо рух уздовж одного кристалічного напрямку, де хвильовий вектор дорівнює \(k\), і виберемо фази наших базових станів таким чином, щоб параметр взаємодії \(P\) був дійсним числом. Наведений нижче гамільтоніан є ілюстративним прикладом:

\[ \boxed{ H(k)= \begin{pmatrix} E_c+\dfrac{\hbar^2 k^2}{2m_0} & Pk\\[6pt] Pk & E_v+\dfrac{\hbar^2 k^2}{2m_0} \end{pmatrix} } \]

Тут, \(P=\hbar\langle u_{c0}|\hat p|u_{v0}\rangle/m_0\) – для обраного напрямку. Зверніть увагу, що два діагональні члени зникають у точці \(k=0\), але збільшуються, коли ми віддаляємося від опорної точки. Отже, модель не описує дві незалежні параболи: стани провідної та валентної зон можуть взаємодіяти.

Щоб отримати допустимі значення енергії, ми розв’язуємо задачу на знаходження власних значень, подібно до того, як ми робили це, коли визначали обмежені стани квантової ями:

\[ H(k)\begin{pmatrix}c_c\\c_v\end{pmatrix} =E(k)\begin{pmatrix}c_c\\c_v\end{pmatrix}. \]

Замість того, щоб знайти єдину енергію, матриця 2x2 створює дві енергетичні гілки. Для цього спрощеного прикладу вони такі:

\[ E_{\pm}(k)= \frac{E_c+E_v}{2}+\frac{\hbar^2k^2}{2m_0} \pm\sqrt{\left(\frac{E_g}{2}\right)^2+P^2k^2}. \]

У точці \(k=0\) обидва результати зводяться до \(E_c\) і \(E_v\), як і очікувалося. У міру збільшення \(k\) взаємодія змінює обидві енергії та змішує основні стани. Ця залежність від квадратного кореня також показує, чому взаємопов’язана зона не обов’язково повинна залишатися ідеально параболічною за межами своєї екстремальної точки.

Практичне заняття важливіше, ніж просто запам’ятовування цієї формули: форма смуги залежить не лише від власної енергії цієї смуги, але й від її взаємодії з сусідніми смугами.

Звідки береться ефективна маса?

У навчальному посібнику про параболічну квантову яму ми використовували ефективну масу як вхідний параметр. Тепер, завдяки наведеному прикладу з двома зонами, ми можемо побачити, звідки вона береться. Дуже близько до значення \(k=0\) результат для зони провідності можна наближено обчислити, розклавши квадратний корінь для малих значень \(k\):

\[ E_+(k)\simeq E_c+ \left(\frac{\hbar^2}{2m_0}+\frac{P^2}{E_g}\right)k^2. \]

Порівняйте це з нашим добре відомим параболічним виразом:

\[ E(k)\simeq E_c+\frac{\hbar^2k^2}{2m_e^{*}}. \]

Обидва мають однакову математичну форму. Порівняння їхніх квадратичних коефіцієнтів дає змогу визначити ефективну масу цієї спрощеної двосмугової моделі:

\[ \boxed{ \frac{1}{m_e^{*}}= \frac{1}{m_0}+\frac{2P^2}{\hbar^2E_g} } \]

Це корисний результат. Ефективна маса не є довільною заміною маси електрона. Вона виникає через те, що стан електрона пов’язаний з іншими електронними зонами всередині кристала. У цьому конкретному випадку менша ширина забороненої зони або сильніша взаємодія призводять до більшої кривизни зони провідності та меншої ефективної маси.

У реального напівпровідника є набагато більше енергетичних зон, ніж лише дві, які ми розглядаємо в нашому прикладі. Якщо ми розглянемо їхню взаємодію як невелике збурення, то відповідна поправка до енергії окремої зони буде такою:

\[ E_n(\mathbf k)\simeq E_n(0) +\frac{\hbar^2k^2}{2m_0} +\frac{\hbar^2}{m_0^2} \sum_{m\ne n} \frac{|\mathbf k\cdot\mathbf p_{nm}|^2} {E_n(0)-E_m(0)}. \]

Кожна сусідня енергетична зона вносить свій внесок відповідно до її імпульсної взаємодії та різниці енергій. Отримана кривина, як правило, залежить від напрямку, тому найбільш повний опис дає тензор ефективної маси, а не одна скалярна маса. Скалярна ефективна маса є доречною, коли симетрія або обраний метод наближення дозволяють розглядати кривину як однакову в обраних напрямках.

Існує одне важливе обмеження: пертурбативний вираз стає ненадійним, коли знаменник малий, тобто коли дві зони знаходяться близько одна до одної за енергією. У такій ситуації краще явно враховувати обидві зони в гамільтоніані, а не приховувати їхню взаємодію всередині єдиної ефективної маси.

Від двох смуг до багатосмугового гамільтоніана

Наступний крок за своєю суттю досить простий. Замість того, щоб починати лише з двох кристалічних станів, ми зберігаємо декілька відповідних станів у нашій основі:

\[ |u_{\mathbf k}\rangle =\sum_{n=1}^{N}c_n(\mathbf k)|u_{n0}\rangle. \]

У результаті отримуємо гамільтоніан розміром \(N\) на \(N\). Кожен діагональний елемент описує опорну зону та її квадратичну поправку; кожен допустимий недіагональний елемент описує взаємодію між двома базисними станами. Розв’язуючи матрицю для обраного хвильового вектора, ми отримуємо всі отримані енергії зон і величину кожного базисного стану у відповідному власному векторі.

Це особливо важливо для валентної зони. Стани важких і легких дірок можуть бути виродженими в центрі зони, і їх змішування неможливо адекватно описати, призначивши кожній із них незалежну постійну масу. Багатозонна модель дозволяє цим станам природно взаємодіяти зі зміною хвильового вектора.

Нам все ще потрібно визначитися, які енергетичні зони слід зберегти. Загальний підхід полягає в тому, щоб залишити ті зони, які найбільше пов’язані з досліджуваними енергіями та фізичними явищами. Взаємодія з більш віддаленими, виключеними зонами може бути приблизно врахована шляхом внесення поправок до параметрів збереженого гамільтоніана. Саме тому різні k·p моделі можуть мати різну кількість енергетичних зон і різні параметри матеріалів.

Наприклад, поширена восьмизонна модель цинкбленду зберігає s-подібний стан провідності та три p-подібні валентні орбітальні стани, кожен з яких має дві спінові компоненти. Детальний вибір базисних функцій, симетрійних членів і елементів матриці пояснюється в нашому розділі базисні функції та гамільтоніани та на сторінках, присвячених конкретним матеріалам.

Навіщо нам потрібні спін і спін-орбітальна взаємодія?

До цього моменту наш простий приклад із двома енергетичними смугами не враховував спін. Однак електрони мають дві складові спіну. Це означає, що кожен квантовий стан повинен бути представлений у вигляді спінових станів «спін вгору» та «спін вниз» (або відповідного еквівалентного спінового базису).

Якби спін не впливав на енергії, це призвело б до повторення проблеми. Реальні напівпровідникові кристали також демонструють спін-орбітальну взаємодію, яка пов’язує спін електрона з його рухом у полі кристалу. У спрощеному описі Паулі додатковий гамільтоніан має такий вигляд:

\[ \hat H_{\mathrm{so}}= \frac{\hbar}{4m_0^2c^2} \boldsymbol{\sigma}\cdot \left(\nabla V_{\mathrm{crystal}}\times\hat{\mathbf p}\right), \]

де \(c\) – швидкість світла, а \(\boldsymbol{\sigma}\) позначає матриці Паулі. Нам не потрібно виводити цей релятивістський член тут; важливим є його фізичний наслідок. Спін і орбітальний момент більше не є незалежними, і певні стани, які в іншому випадку мали б однакову енергію, розділяються.

У добре відомій валентній зоні структури цинкбленду спін-орбітальна взаємодія призводить до утворення важких, легких та розділених підзон, причому енергія розділеної підзони зазвичай позначається як \(\Delta_{\mathrm{so}}\). Це не просто три непов’язані параболи: їхні властивості та взаємодія залежать від хвильового вектора та основної кристалічної симетрії.

Це пояснює назву «восьмизонна k·p модель»: чотири орбітальні базові стани (один провідний і три валентні), кожен з яких має дві спінові складові, утворюють восьмикомпонентну базу. Фактичний гамільтоніан включає дозволені значення імпульсу та спін-орбітальні взаємодії між ними. Повний розрахунок для квантової ями з використанням цієї моделі продемонстровано в навчальний посібник з восьмизонної цинкової структури k·p для GaAs/AlGaAs.

Як це перетворюється на обчислення, що використовує модель квантової ями?

До цього моменту ми описували об’ємний кристал, у якому енергія електронів змінюється залежно від тривимірного хвильового вектора. Щоб описати квантову яму, нам тепер потрібно повернутися до геометрії, представленої на попередній сторінці навчального матеріалу.

Припустимо, що напівпровідникові шари розташовані паралельно площині \(x\)-\(y\) і утворюють стовпчик вздовж осі \(z\). Вектор хвилі в площині, \(k_{\parallel}\), залишається корисним параметром, оскільки структура приблизно однорідна вздовж шарів. Однак, електрон обмежений у напрямку \(z\), тому нам потрібно розв’язати рівняння для визначення його просторової зміни в межах області.

У стандартному підході, що базується на використанні амплітудної функції, вектор хвилі, який визначає напрямок росту, у гамільтоніані k·p замінюється диференціальним оператором:

\[ \boxed{k_z\;\longrightarrow\;-i\frac{d}{dz}} \]

Межі зон і інші відповідні параметри матеріалу також можуть змінюватися залежно від \(z\), що відображає різні напівпровідникові шари. Проблема полягає не просто в обчисленні звичайної скінченної матриці для одного хвильового вектора: це набір пов’язаних диференціальних рівнянь для амплітудних функцій у всьому пристрої. Коли параметри змінюються на межах, порядок застосування операторів і умови узгодження мають бути обрані послідовно, щоб гамільтоніан залишався фізично значущим і ермітовим.

Найпростіша параболічна модель розв’язує одне скалярне рівняння Шредінгера для електрона або дірки. Натомість багатосмугова модель k·p розв’язує систему рівнянь для кількох взаємопов’язаних компонент амплітуди. Отримані власні значення визначають обмежені енергетичні рівні, а їхні власні вектори показують, як збережені характеристики смуги впливають на кожен стан.

Повторюючи обчислення для різних значень \(k_{\parallel}\), ми формуємо структуру субзон квантової ями. На відміну від набору незалежних парабол, ці субзони можуть взаємодіяти, мати непараболічну форму та змінювати свої властивості зі збільшенням хвильового вектора в площині.

Що саме ми можемо обчислити завдяки додатковим знанням з фізики?

Перевага багатосмугових обчислень полягає не лише в тому, що вони дозволяють отримати складніші діаграми зон. Вони дають нам більш реалістичний опис електронних станів, які беруть участь у транспортних і оптичних процесах.

У нашому навчальному посібнику з оптичних переходів ми розглянули ймовірність переходу електрона між заповненим і незаповненим станом. Інтенсивність переходу залежить від матричного елемента між початковою та кінцевою хвильовими функціями. Багатокомпонентний власний вектор містить інформацію про змішаний орбітальний і спіновий характер кожного стану, що є важливим під час обчислення напрямку та інтенсивності оптичних переходів.

Це особливо корисно під час вивчення переходів, що залежать від поляризації, оптичного підсилення, спонтанного випромінювання та систем, у яких спостерігається значне змішування валентних зон. Детальні оптичні матричні елементи залежать від обраного гамільтоніана та бази, але їх обчислення природно випливає після отримання пов’язаних власних станів.

k·p моделі в OghmaNano

OghmaNano пропонує декілька моделей для різних напівпровідникових структур. Усі вони базуються на загальній філософії k·p, розробленій тут, але мають різні електронні стани та використовують матеріальні параметри, що відповідають кожній кристалічній системі.

• Восьмизонна модель цинкбленду k·p включає провідні та p-подібні валентні стани, а також спін і спін-орбітальну взаємодію.
• Десятисмугова модель цинкбленду k·p зберігає розширену електронну основу для розв’язання задач, які потребують додаткових взаємопов’язаних станів.
• Восьмизонна модель Вурцита k·p використовує симетрію та параметри зон, які відповідають кристалам Вурцита.

Цей. наступний розділ описує фактичні базові стани та пояснює, як вони організовані у вигляді гамільтонових матриць, які використовуються цими розв’язувачами.