1 МІНІСТЕРСТВО ОСВІТИ I НАУКИ УКРАЇНИ Харківський національний університет імені В. Н. Каразіна Факультет радіофізики, біомедичної електроніки та комп'ютерних систем Кафедра фізичної і біомедичної електроніки і комп'ютерних інформаційних технологій ЗАТВЕРДЖУЮ В.о. завідувача кафедри _________Сергій_БЕРДНИК_ підпис ініціали, прізвище “____” ____________2019 року Кваліфікаційна робота магістра на тему: ВАРІЗОННІ ДІОДИ ДЛЯ ГЕНЕРАЦІЇ КОЛИВАНЬ В ТЕРАГЕРЦОВОМУ ДІАПАЗОНІ Виконав: студент ІІ курсу магістратури, групи РE-61 спеціальності 176 Мікро- та наносистемна техніка освітньо-професійна програма «Фізична та біомедична електроніка» Гребенюк Р. Ю. Керівник кандидат фіз.-мат. наук, доцент Боцула О.В. 2023 рік 2 РЕФЕРАТ Дипломна робота містить 54 сторінок, 20 ілюстрації, 1 таблиця, 49 бібліографічних джерел. В роботі проведено дослідження електромагнітних коливань в терагерцового діапазону варізонними діодами на основі InPAs з використанням моделювання на основі багаточасткового методу Монте – Карло з урахуванням ударної іонізації. Показано, що ВАХ коротких варизонних діодів не містить області з негативною диференціальною провідністю, однак вони демонструють нестійкості струму типу зарядженого шару та генерацію коливань. Дослідження показало, що максимальна ефективність генерації становить приблизно 10%, спостерігається в діодах довжиною 1280 нм на частоті 100 ГГц. У більш коротких діодах ефективність знижується до 3,9% і 2,0% у діодах довжиною 640 і 500 нм відповідно. Гранична частота генерації була близько 400 ГГц в діодах довжиною 500 нм. КЛЮЧОВІ СЛОВА: терагерцовий діапазон, діод, ударна іонізація, варізонний шар, ефективність генерації 3 ABSTRACT This work contains 54 pages, 20 illustrations, 1 table, 49 bibliographic sources. А study of electromagnetic oscillations in the terahertz range by graded-band diodes based on InPAs has been carried out using the Ensemble Monte Carlo Technique with consideration of impact ionization. It was shown that the I-V characteristic of short graded-band diodes does not contain areas with negative differential conductivity, but they demonstrate charged layer current instabilities and oscillation generation in resonant circuits. The study revealed that the maximum generation efficiency is approximately 10%, observed in diodes with a length of 1280 nm at a frequency of 100 GHz. In shorter diodes, the efficiency decreases to 3.9% and 2.0% in diodes with lengths of 640 and 500 nm, respectively. KEYWORDS: terahertz range, diode, impact ionization, graded-gap layer, generating efficiency, 4 ЗМІСТ 1. ОГЛЯД ЛІТЕРАТУРИ .................................................................................... 7 1.1 Ударна іонізація в налвисокочастотних приладах як фактор підвищення їх частотних можливостей ............................................................ 7 1.2 Ударна іонізація в діодах на основі сполуки Gal - x InxAsyP1 - y 1.3 Математичне моделювання напівпровідникових приладів ............... 21 1.3.1 Дрейфово-дифузійна модель ......................................................... 24 1.3.2 Гідродинамічні моделі 26 1.3.3 Метод Монте – Карло ............................................................................ 27 2. МОДЕЛЮВАННЯ РОБОТИ ВАРІЗОННОГО ДІОДА НА ОСНОВІ InPAs ....................................................................................................................... 30 2.1 Модель діода 32 2.1.2 Вразуовання розсіяння носіїв заряду ................................................... 32 2.1.3 Алгоритм методу Монте-Карло ............................................................. 36 2.2. Результати моделювання 39 2.2.1 Постановка зададачі 39 2.2.2 Характеристки діодів на постійному струмі 42 2.2.3 Частотні характеристики діодів 44 ВИСНОВКИ 49 ПЕРЕЛІК ПОСИЛАНЬ 50 5 ВСТУП Зі збільшенням поширення ТГц технології одним із важливих завдань є розробка компактних джерел випромінювання, придатних для інтеграції в компактні терагерцові системи [1]. Важливу роль у цьому напрямку відіграють резонансні тунельні діоди (RTD) з частотою зрізу понад 2 ТГц [2,3]. Однак недоліком RTD є їх відносно низька вихідна потужність. Останнім часом існує оптимізм щодо використання квантово-каскадних лазерів (ККЛ) [4,5], які є перспективними джерелами для різних додатків ТГц. Однак основним недоліком ККЛ є необхідність додаткового охолодження. Максимальна робоча температура ТГц ККЛ становить 250 К в імпульсному режимі [5]. Потенційним підходом є покращення характеристик звичайних пристроїв, таких як діоди Ганна або діоди з часом лавинного проходження (ATTD). Діоди Ганна демонструють значне зниження потужності змінного струму на частоті, корельованій із зворотним міждолинним часом перенесення електронів. Наприклад, використання бінарних сполук, таких як GaN, InN і GaN, а також потрійних III-нітридних сполук може давати пристрої з покращеними частотними характеристиками. Однак експериментальна реалізація діодів Ганна на основі нітриду галію залишалася невирішеною протягом тривалого періоду. У [6] запропоновано діод Ганна на основі GaN довжиною 0,6 мкм, де отримано коливання на основній частоті в діапазоні 0,3 – 0,4 ТГц. На сьогоднішній день значного покращення властивостей пристроїв можна досягти за рахунок використання гетероструктур, неузгоджених матеріалів та сортованих напівпровідників [7]. Застосування градієнтних напівпровідників у пристроях з перенесеними електронами покращує умови для транспорту електронів як на контактах, так і в області проходження діода [8, 9]. Діоди на основі градієнтних напівпровідників можуть перевершувати однорідні діоди як за ефективністю генерації, так і за частотними 6 властивостями. Особливий інтерес представляють випадки, пов’язані з вузькозонними напівпровідниками з малими значеннями енергетичної щілини, такими як потрійні III-V сплави, такі як GaInAs, AlInN, InPAs, InGaN. Тут InN та InAs демонструють максимальні швидкості та мінімальне порогове електричне поле для переносу електронів із нижньої долини зони провідності у верхні сателітні [10,11]. Використання вузькозонного напівпровідника на анодному контакті збільшує ширину забороненої зони між нееквівалентними долинами. Крім того, ці напівпровідники мають низьку порогову енергію для ударної іонізації (II), що робить неможливим виготовлення діодів на основі InN або InAs з цієї причини. Однак у діодах із градієнтними шарами обмежена ударна іонізація може бути чинником для покращення частотних властивостей. Це явище було продемонстровано на діоді з градиентним шаром InGaAs [12,13]. У варизонному напівпровіднику n-типу зона провідності плоска. Отже, при прикладеному зміщенні ударна іонізація ініціюється лише електронами. Крім того, при помірних прикладених напругах зміщення електрони та дірки можуть рухатися в одному напрямку до анода. При вищих прикладених напругах деякі дірки переміщуються до катода, але ударна іонізація, ініційована дірками, не відбувається через збільшення енергетичного зазору та малу силу, що діє на дірку. Отже, на аноді відбувається релаксація енергії електронів. Показано, що в діоді довжиною 500 нм частота генерації досягає 320 ГГц [13]. Іншою можливою альтернативою GaAs як основного матеріалу є використання InP і відповідна заміна сортованого InGaAs сортованим InPA. InP має вищу швидкість електронів у порівнянні з GaAs і демонструє більший пробій за тієї самої напруги зміщення [14]. Метою даної дипломної роботи є дослідження генерації електромагнітних коливань у довгохвильовій частині терагерцового діапазону за допомогою варизоннлшл діода на основі системи матеріалів InP/InPAs. 7 1. ОГЛЯД ЛІТЕРАТУРИ 1.1 Ударна іонізація в налвисокочастотних приладах як фактор підвищення їх частотних можливостей Найважливішим елементом модемної ТГц технології є генератор когерентного ТГц випромінювання. Напівпровідникові джерела представляють великий інтерес завдяки компактним розмірам і можливості інтеграції в один кристал. Найпоширенішим джерелом є діоди, що працюють на ефекті між долинного переносу електронів, діоди Ганна. Підвищення частоти генерації цих пристроїв може бути досягнуто лише шляхом зменшення розмірів пристрою, що в основному призводить до появи дуже сильного електричного поля всередині пристрою під час роботи. У таких умовах ударна іонізація може призвести до неконтрольованого зростання струму, що може призвести до виходу нз ладу діода. З цієї точки зору великий інтерес представляє верхня межа частоти стаціонарного ефекту Гімна в пристроях, що використвувються для генерації, діодах та транзисторах. Ефект Гімна в умовах ударної іонізації в GaAs досліджувався, наприклад, в MESFET, де проводилося моделювання процесів переносу заряду методом частинок, Монте-Карло (MC'P) [15]. Було показано, що в таких структурах можна отримати коливання струму стоку на чатості 125 ГГц при довжині затвора 400 нм. Також це моделювання з використання методу Монте-Карло показало виникненя коливання струму стоку на частоті 50 ГГц в умовах ударної іонізації у польовму трнзистрів з високою рухливістю електронів AlGaAs/InGaAs [16]. В роботі [14] цей процес було розглянуто в МОП-транзисторах на основі InP за умов ударної іонізації та шляхом правильного підбору параметрів структури, напруги зміщення, який 8 дав змогу отримати найвищі частоти коливань струму стоку, які виникали внаслідом ефекту Ганна. Результати цих розрахунків для структури 50-L-50 нм n+nn+- на основі InP каналу шириною 200 нм показані на рис. 1.1 Рис. 1.1 Залежність струму витоку від напруги на витоку прирізним напругах на затворі[14] Рис. 1.2 Залежність струму від часу[14] Довгий затвор Lz центруваася в n-області на відстані 50 нм від каналу. Легування n = 5×10 16 см −3 і n+ =4×10 18 см −3 . Моделювання було виконано шляхом одночасного розв’язання узгоджених рівнянь Больцмана та двовимірного рівнянь Пуассона [17]. Зона провідності InP описувалася J, A·10 9 А/м 2 U, В J, A·10 9 А/м 2 t, пс 9 трьома долинами провідності (Γ-L-X) [17], валентна - двома валентними (зона важких дірок та спін-орбітально відщепленої) зонами[18]. Також модель включала ударну іонізацію, час життя носіїв і Оже - рекомбінацію [19]. Температура решітки в усіх випадках становила 300 К. Результати, наведні на рис. 1.1. та рис.1.2 відповідвають для структурі з L =0,7 мкм та Lz=0,5 мкм. Стаціонарні ганнівські rколивання струму стоку на частоті 200 ГГц виникають при достатньо високому (приблизно 3 В) прямому зміщенні затвора, де зміщення стоку UC становить 1,8 В (рис. 1.2). Ударна іонізація у цьому випадку відсутня. Вона виникає при UC понад 2,0 В. За рахунок ударної іонізації амплітуда коливань стаціонарного струму стоку Ганна та середнє значення струму зростають із збільшенням UC до 3,1 В (рис.1.2). Амплітуда струму при UC = 3,1 В перевищує величину, яка була розрахована без урахування ударної іонізації. Амплітуда струму зростає за рахунок синхронізації інтенсивності ударної іонізації з коливаннями струму (рис.1.3). Рис. 1.3 Залежнісь густини струму від часу протягом чотирьох періодів коливань струму стоку (суцільна лінія) та числа актів ударної іонізації (пунктирні лінії). J, A·10 9 А/м 2 t, пс 10 На рис 1.3 пронумеровані точки від 1 до 4 позначають моменти часу, в які обчислюється концентрація та енергія електронно-діркових пар. Напруги на стоці і затворі відпвідно UC = 3,1 В і UЗ = 3 В Максимальний струм і число актів іонізації досягається, коли домен Ганна входить у контакт стоку. Утворення та дрейф домену Ганна відбуваються при дуже низькому рівні ударної іонізації. Згенеровані дірки через високі ефективні маси не беруть участі в коливаннях струму стоку Ганна. Вони служать додатковим легуванням каналу. Було показано, що максимальна стаціонарна частота коливань струму стоку Ганна складає 500 ГГц і вдповідає ПТ на основі InP з L = 0,4 мкм і LЗ = 0,3 мк. Мінімальна напруга затвора для коливань Ганна струму стоку становить 8 В. Коливання починаються при UC = 2 В. При UC >2,8 В відбувається неконтрольоване зростання струму. Ефект виникнення коливань найбільше проявляється при позитивних напругах на затворі. При позитивному зміщенні електрони з високолегованих областей витоку та стоку промують в бік затвору, залишаючи нижчу концентрацію електронів у областях за межами затвора. Концентрація електронів в області між витоком і затвором залишається стабільною навіть під час виникнення коливань (рис.1.4). дірки електрони x, мкм n,p,·10 17 ·м 3 11 Рис. 1.4 Розподіли концентрацій електронів і дірок у різні моменти часу протягом періода коливань, які відпровідають рис.1. 3. Концентрація дірок показана в нижній частині зображення. Це призводить до утворення високого поля в цій області. Поле запускає електрони високої енергії в канал під затвором, що сприяє формуванню домену Ганна. В роботі[20] дослідження, аналогічні попереднім, проведені для транзисторної уніполярної структури, яка являла собою 50-320-50 нм n+nn+ GaAs канал з довжиною 200 nm. Затвор довжиною 200 нм позміщувался на n-обьласті на відстані 80 нм від n+- витоку. Ріень легування в каналі вибраний рівним 5·10 22 м -3 та для n+ 10 25 м -3 . Ефективна відстань між завтором і каналом вибиралася рівною 100 нм. Модель для опису зон вибьрана як і попередньому випадку. Крім того, в модель включили ударну іонізацію, кА описувалася формулою Келдиша та враховувалася Оже – рекомбінація з точним коефіцієнтом для електронів і кульок C = 10 -28 см 6 с -1 [21]. Мінімальна температура в усіх випадках є 300 K. Кількість частинок 400,000, а проміжок часу ї 0.5 фс. Лавинний пробій у структурі відбувається при напрузі стоку понад 4,5 В, а процеси ударної іонізації починаються понад 2,6 В. Струм стоку починає зростати при напрузі стоку 3,8 В, коли інтенсивність ударної іонізації є достатньо високою. Тому за співвідношенням струм-напруга важко ідентифікувати початок процесів ударної іонізації. Струми стоку як функції часу при різних зміщеннях стоку показані на рис. 1.5 Знизу вгору: зміщення діоду UC = 2, 3 і 4 В. Щоб уникнути перекриття, криві зсуваються на 0, 0,5 і 1,5·10 9 A·M 2 . На рис.1.3 можна побачити низьку амплітуду коливань струму стоку при UC = 3 В, коли ударна іонізація незначна. Коливання великої амплітуди виникають при UC = 4 В. У цьому випадку процеси ударної іонізації є інтенсивними. Висока амплітуда коливань струму стоку виникає за рахунок ударної іонізації та синхронізації руху домену Ганна. Максимальний 12 струм і кількість подій ударної іонізації досягаються, коли домен Ганна входить у контакт стоку. Спектральні щільності струму стоку при різних зміщеннях стоку наведені на рис.1.5 Для більш детального аналізу процесів у каналі приведені розподіли, показані на рис. 1.7 Рис. 1.5 Залежнісь густини струму від часу в залежності від величини напруги витік - стік при напру Очевидний пік основного коливання струму стоку на частоті 325 ГГц (рис.1.6). Рис. 1.6 Спекральна густина струму як функція частоти при різних значеннях напруги витік - стік [20] J, в.о. f, ГГц J, A·10 9 А/м 2 t, пс 13 Крім того, при зміщенні стоку 4 В можна побачити піки, що відповідають другій та третій гармонікам. У цьому випадку з'являється посилений низькочастотний шум, який походить від діркової складової струму. Через велику ефективну масу та низьку швидкість у каналі (рис. 1.7, б) дірки не впливають на основну пікову частоту, як це можна побачити порівнюючи піки при 3 В та 4 В на рис. 1.6. На рис. 1.7, a) концентрація електронів у каналі під затвором є сильною посиленими порівняно з рівнем легування, а біля витоку утворюється область зниженої концентрації. ЕЛЕКТРОНИ ДІРКИ x, нм x, нм а) x, нм ЕЛЕКТРОНИ ДІРКИ б) в) n,p,·10 22 ·м 3 v,·10 5 ·м/c E, кВ/см 14 Рис. 1.7. Розподіли концентрації електронів і дірок (а), швидкості (б) і електричного поля (в) при зміщенні стоку 4 В і UЗ = 25 В У діодах Ганна виїмка виконується нижчим легуванням області поблизу джерела. Концентрація дірок у каналі менша на два порядки порівняно з концентрацією електронів і експоненціально розпадається в контактах витоку та стоку через рекомбінацію Оже. Отже, дірками в динаміці ефекту Ганна можна знехтувати. Рисунок 1.7 показує рух шару накопичення електронів у каналі. Видно, що накопичувальний шар з’являється ближче до правого боку виїмки і переміщується до контакту витоку. Еволюція енергії електронів під час генерації Ганна продемонстрована на рис. 1.8 Рис. 1.8 Просторові розподіли енергії електронів у чотори різні моменти часу з інтервалом близько часу 0,6 пс. Напруги на транзисторі UC =4 В і UЗ = 25 В. [20] Стабільна підвищена енергія в області надрізу під час генерації очевидна. Це типово для роботи діода Ганна з вирізом. Максимальна енергія та інтенсивність ударної іонізації досягаються при вході накопичувального шару в контакт витоку. Це є причиною підвищення амплітуди струму стоку в умовах ударної іонізації. Еволюція профілю електричного поля під час генерації Ганна показана на рис. 1.9 E, eB x, нм 15 Fig. 1.9 Просторові розподіли напруженості електричного поля у чотори різні моменти часу з інтервалом близько часу 0,6 пс. Напруги на транзисторі UC =4 В і UЗ = 25 В. Отже, ефект Ганна в польових транзисторах є найбільш ефективним при позитивних напругах на затворі. При позитивному зміщенні електрони з високолегованих контактах витоку та стоку захоплюються затвором, що приводить до зменщення концентрація електронів у прелеглих до нього областях. Концентрація електронів в області між витоком і лівим краєм затвора залишається стабільною навіть під час генерації. Це призводить до формування стабільного високого поля в цій області (рис. 1.9). Поле створює високоенергетичні електрони в канал під затвором, що спияє формуванню домену Ганна або зарядженого шару. Результатом є виникнення коливань високої частоти, для польовому транзисторі на основі GaAs вона може досягати величини 325 ГГц. 1.2 Ударна іонізація в діодах на основі сполуки Gal - x InxAsyP1 - y Результати експериментального дослідження ударної іонізації в діодах на основі сполуки Gal - xInxAsyP1 - y , яка мала орієнтацію (100), наведено в роботі[22]. Характеристки отримані з вимірювань фотопомноження, які E, кВ/см x, нм 16 проводилися на п’яти видах діодів Gal-xInxAsyP1-y/InP. Ці залежності аналізувалися з використанням теорії Барафа щодо зв’язку між коефіцієнтом іонізації та електричним полем. Було проведено дослідження також температурних залежностей напруги пробою діодів InP і GalnAsP. У цьому дослідженні використовуються діоди меза-типу /InP, вміст миш’яку яких становить 0, 0,28 , 0,39 , 0,70 і 1,0. Було виготовлено три типи діодів (типи від А до С). Вони схематично показані на рис. 1 [(a) тип A, (b) тип B, (c) тип C]. Рис. 1.10 Схематичний перереріз Gal-xInxAsyP1-y/InP - діодів Тип А використовувався для діода I, тип B для діодів II - IV і тип C для діода V відповідно. Тип А використовувався для діода I, тип B для діодів II - IV і тип C для діода V відповідно. Структури діодів і параметри пристрою наведені в таблиці 1. Діод 1 використовує шар зупинки травлення n-GaInAsP для вибіркового травлення підкладки - InP. Діоди 1-4 були вирощені методом рідкофазної епітаксії, в якій використовувався звичайний метод горизонтального ковзання. Лише діод 5 був вирощений методом епітаксії з А Б С 17 парової фази (VPE), оскільки в його допомогою можна легко вирощувати n- шари, які необхідні особливо в GaInAs/InP діодах, щоб мінімізувати ефект міжзонного тунелювання. Таблиця 1 № діода 1 2 3 4 5 Матеріал діода InP GalnAsP (y-0.28) GalnAsP (y-0.39) GalnAsP (y=0.70) GalnAs Структура A Б Б Б В Товщина активної частини, мкм 2,1 3,0 1,7 1,7 2,3 Концентрація в активній частині n·10 14 , см -3 300 100 170 200 10 Товщина n-InP буферного шару, мкм — 6,0 4,8 5,0 4,8 Концентрація в n-InP, n·10 14 , см -3 100 170 200 5 У всіх діодах шляхом дифузії Cd утворювалися -шари глибиною 1,3- 1,5 пм. Суміш Br:HBr:H20 у співвідношенні 1:17:34 використовували як мезатравник. Омічні контакти до - та - шарів створювали шляхом легування напиленого AuZn та AuGe відповідно. Діаметри світлочутливої зони та ділянки діода становлять відповідно 100 і 200 пм. Вікно діаметром 80 мкм в електроді з боку підкладки було сформовано для того, щоб опромінювати світло із задньої частини діода. Діоди 1 та 2 показали характеристики лавинного пробою, а темновий струм при 90% напруги пробою становив лише 10 нА. У діоді 3I тунельний струм виявився близьким до пробою, але він був менше приблизно 20 пА, оскільки ширина забороненої зони кристала GaInAsP була 18 відносно широкою. У діодах 4 і 5 пробій стабілітрона був домінуючим при високих напругах зворотного зміщення через вузьку заборонену зону. Залежності від напруги зворотного зсуву коефіцієнтів розмноження, ініційованих електроном і діркою, , і , були виміряні шляхом випромінювання світла з довжинами хвилі і з поверхні і задньої частини кожного діода відповідно. Зазначені довжини хвиль і відбиралися монохромометром, який освітлювався вольфрамовою лампою. Типові характеристики фотомноження діодів показані на рис. 1.11 Рис. 1.11 Характеристики фотопомноження діодів Використана довжина хвилі становить 0,83 пм для всіх діодів, а довжини хвилі становлять 0,63 , 0,97 , 0,97 , 1,15 і 1,25 пм для діодів I – V відповідно. 19 Встановлено залежності коефіцієнтів помноження для електронів та доірок від напруженості електричного поля для кристалів п’ятох типів сплаву Gal-xInxAsyP1-y. , що показано на рис. 1.12. Рис. 1.12 Залежність коефіцієнтів ударної іонізації електронів та дірок від величини зворотної напруженості електричного поля Ці графіки чітко вказуть на залежність коефіцієнтів ударної іонізації для електронів  та дірок  від молярного вмісту миш'яку. В діапазоні 0 1. Припумкається,що ця зміна може бути обумовлена двома причинами: одна - це ефект розсіювання, пов’язаного з розсіянням на сплаві, яке впливає на енергетичні дірки більше, ніж на енергетичні електрони в кристалах 21 GaInAsP. Друга - це ефект набагато меншого ефективного часу розсіювання для енергетичних електронів при великому вмісті миш'яку. 1.3 Математичне моделювання напівпровідникових приладів Оскільки розміри напівпровідникових елементів зменшуються до режиму нанометрового масштабу, навіть поведінка звичайних пристроїв стає дедалі складнішою, оскільки виникають нові фізичні явища при малих розмірах і досягаються обмеження у властивостях матеріалу [23]. На додаток до проблеми, пов’язані з розумінням фактичної роботи надмалих пристроїв, зменшені розміри функцій вимагають більш складних і трудомістких процесів виробництва. Цей факт означає, що чистий підхід до оптимізації пристрою методом проб і помилок стане неможливим, оскільки він забирає надто багато часу та занадто дорогий. Оскільки комп’ютери є значно дешевшим ресурсом, симуляція стає незамінним інструментом для розробників пристроїв. Крім того, пропонуючи можливість тестувати гіпотетичні пристрої, яких ще немає (або не міг) ще бути виготовлений, симуляція пропонує унікальне розуміння поведінки пристрою, дозволяючи спостерігати за явищами, які неможливо виміряти на реальних пристроях. Обчислювальна електроніка [24,25] у цьому контексті відноситься до фізичного моделювання напівпровідникових пристроїв з точки зору транспортування заряду та відповідної електричної поведінки. Це пов’язано, але зазвичай окремо від моделювання процесу, яке має справу з різними фізичними процесами, такими як ріст матеріалу, окислення, дифузія домішок, травлення та осадження металу, властиві виготовленню пристроїв [26], що веде до інтегральних схем. Симуляцію пристроїв можна розглядати як один із компонентів технології автоматизованого проектування (TCAD), яка забезпечує основу для моделювання пристроїв, яке має справу з 22 компактними поведінковими моделями для пристроїв і підсхем, що мають відношення до моделювання схем у комерційних пакетах. такі як SPICE [27]. Головна мета обчислювальної електроніки полягає в тому, щоб створити інструменти моделювання з необхідним рівнем складності для охоплення суттєвої фізики, водночас мінімізуючи обчислювальне навантаження, щоб результати могли бути отримані в розумні часові рамки. Є дві основні задачі, які необхідно розв’язувати узгоджено одна з одною, які пов’язані з описом процесів переносу, що описують потоки заряду, і знаходженням поля, що керує потоком заряду. Обидві проблеми тісно пов’язані одна з одною, тому їх потрібно вирішувати одночасно. Великі поля виникають із зовнішніх джерел, а також від щільності заряду та струму, які діють як джерела змінних у часі електричних і магнітних полів, отриманих у результаті розв’язання рівнянь Максвелла. За відповідних умов діють лише квазістатичні електричні поля, що виникають внаслідок рішення рівняння Пуассона. Поля, у свою чергу, є рушійними силами для транспортування заряду, Вони визначаються для різних рівнів наближення в ієрархічній структурі, починаючи від компактного моделювання у верхній частині до точного квантово-механічного опису в нижній частині. На самому початку напівпровідникової технології характеристики електричних пристроїв можна було оцінити за допомогою простих аналітичних моделей (поступове наближення каналу для МОН-транзисторів), спираючись на формалізм дрейф-дифузії (ДД). Необхідно було зробити різні наближення, щоб отримати рішення закритої форми, але отримані моделі відобразили основні характеристики пристроїв [28]. Ці наближення включають спрощені профілі легування та геометрію пристроїв. У зв’язку з постійним удосконаленням і вдосконаленням технології ці наближення втратили свою основу, тому потрібен був більш точний опис. Цієї мети можна досягти шляхом чисельного вирішення рівнянь. Чисельне моделювання переносу носіїв у напівпровідникових пристроях сходить до 23 відомої роботи Шарфеттера та Гуммела [29], які запропонували надійну дискретизацію рівнянь дрейфу – дифузії, яка використовується досі. Марківський процес у часі, метод функцій Грінса дозволяє одночасно розглядати кореляції в просторі та часі, обидва з яких, як очікується, будуть важливими в нанорозмірних пристроях. Однак виникають труднощі в розумінні різних термінів у результаті рівняння та величезне обчислювальне навантаження, необхідне для його фактичного впровадження, роблять корисним розуміння квантових ефектів у реальних пристроях обмежених значень. Наприклад, єдиним успішним комерційним використанням підходу функції Гріна є симулятор NEMO (Nano-Electronics Modeling) [30], який фактично є еталоном, і в першу чергу застосовується до резонансних тунельних діодів. Зі сказаного вище випливає, що, всупереч останнім технологічним досягненням, сучасний рівень техніки в моделюванні пристроїв наразі не дає можливості розглядати ці нові проблеми у масштабуванні розмірів пристроїв від звичайних до пристроїв квантового масштабу. Для кремнієвих пристроїв з активними областями менше 0,2 мкм в діаметрі описи макроскопічного транспорту на основі моделей дрейфу-дифузії не можна застосовувати. Як уже зазначалося, навіть стандартні гідродинамічні моделі зазвичай не забезпечують достатньо точний опис, оскільки вони нехтують значними внесками від хвоста функції розподілу фазового простору в областях каналу [31,32]. У рамках вимоги самоузгодженого вирішення проблеми пов’язаного транспортного поля в цій новій області фізики пристроїв є кілька обчислювальних проблем, які обмежують цю здатність. Одна з них полягає в необхідності розв’язувати, як транспортне, так і рівняння Пуассона в усій 3D області пристрою (і за її межами, якщо врахувати ефекти випромінювання). Як наслідок, високоефективні алгоритми, націлені на обчислювальні платформи високого класу (швидше за все, у багатопроцесорних середовищах) необхідні для повного вирішення навіть відповідних польових 24 завдань. Відповідний рівень наближення, необхідний для охоплення належної фізики нерівноважного транспорту, що стосується майбутньої моделі пристрою, є ще більш складною проблемою як з точки зору обчислень, так і з точки зору фундаментальної фізики 1.3.1 Дрейфово-дифузійна модель Для традиційного моделювання напівпровідникових пристроїв переважна модель відповідає розв’язкам так званих рівнянь дрейфу-дифузії, які є «локальними» з точки зору рушійних сил (електричних полів і просторових градієнтів густини носіїв), тобто струм у певній точці простору залежить лише від миттєвих електричних полів і градієнта концентрації в цій точці. Повна дрейфово-дифузійна модель базується на такому наборі рівнянь беперервності є законами збереження носіїв. Числова схема, яка розв’язує рівняння неперервності, повинна 1) Зберігати загальну кількість частинок всередині моделюваного пристрою. 2) Зважати на локальну позитивно визначену природу щільності носія. Від’ємна щільність є нефізичною. 3) Дотримуватися монотонності рішення (тобто воно не повинно вводити помилкових просторових коливань). Консервативні схеми зазвичай досягаються шляхом поділу обчислювальної області на ділянки (комірки), що оточують точки сітки. Потім струми визначаються на межах цих елементів, таким чином забезпечуючи збереження (струм, що виходить з одного елемента сторона точно дорівнює струму, що надходить у сусідній елемент через спільну сторону). За відсутності умов генерації-рекомбінації єдиний внесок у загальний струм пристрою виникає від контактів. Треба пам’ятати, що 25 оскільки електрони мають негативний заряд, потік частинок протилежний потоку струму. Рівняння дискретизуються, наприклад, використовуючи кінцеві різниці, існують обмеження щодо вибору розміру сітки та часового кроку [33]: 1) Розмір сітки Δx обмежений довжиною Дебая. 2) Крок за часом обмежений часом діелектричної релаксації. Розмір сітки має бути меншим за довжину Дебая, коли потрібно розв’язувати зміни заряду в просторі. Простим прикладом є перерозподіл носіїв на межі розділу між двома областями з різними рівнями легування. Носії дифундують у нижню леговану область, утворюючи надлишковий розподіл носіїв, який у стані рівноваги спадає в просторі до об’ємної концентрації з приблизно експоненціальною поведінкою. Константа просторового розпаду є довжиною Дебая. Рівняння дрейф-дифузії напівпровідника складають пов’язану нелінійну систему. Загалом неможливо отримати рішення безпосередньо за один крок, але потрібен метод нелінійної ітерації. Двома найбільш популярними методами розв’язання дискретизованих рівнянь є ітераційний метод Гуммеля [34] і метод Ньютона [35]. Дуже важко визначити оптимальну стратегію для вирішення, оскільки це буде залежати від ряду деталей, пов'язаних з конкретним пристроєм, що досліджується. Нарешті, дискретизація рівнянь нерозривності у формі збереження вимагає визначення струмів у середніх точках ліній сітки, що з’єднують сусідні вузли сітки. Оскільки розв’язки доступні лише у вузлах сітки, то для визначення струмів потрібні схеми інтерполяції. Підхід Шарфеттера і Гуммеля [29] забезпечив оптимальне вирішення цієї проблеми, хоча математичні властивості запропонованої схеми були повністю визнані набагато пізніше. У спільноті обчислювальної електроніки необхідність гідродинамічної (ГД) транспортної моделі зазвичай перевіряється шляхом порівняння результатів моделювання для моделювання ГД і ДД. Незважаючи на 26 очевидний факт, що залежно від набору рівнянь враховуються різні основні фізичні ефекти, вплив на моделі для фізичних параметрів більш тонкий. Основною причиною цього є те, що у випадку моделі HD інформація про середню енергію носія доступна у формі температури носія. Багато параметрів залежать від цієї середньої енергії носія, наприклад, рухливість і час релаксації енергії. У випадку моделі ДД передбачається, що температури носіїв знаходяться в рівновазі з температурою решітки, тобто TC =TL, отже, усі енергозалежні параметри повинні бути змодельовані іншим способом. 1. 3.2 Гідродинамічні моделі Перші три рівняння балансу, отримані шляхом взяття моментів переносу Больцмана. Рівняння (КРБ), набуває вигляду: 1 n n n S t q      ( 1.1) 2 0 2 1 ( )in w i i i m W q W nq F E E W W t m x m                ( 1.2) 22 1n in i n i i m J q W nq E t m x m            ( 1.3) Рівняння балансу для густини носія вводить густину струму носія, яке балансове рівняння вводить густину кінетичної енергії. Рівняння балансу для густини кінетичної енергії, з іншого боку, вводить потік енергії. Таким чином, в ієрархії рівнянь балансу з’являється нова змінна, і набір нескінченних рівнянь балансу фактично є розв’язком БТЕ. Імпульс і швидкості релаксації енергії, які відображаються в рівнянні. – усереднені по ансамблю величини. Для простих механізмів розсіювання можна використовувати дрейфовe максвеллівська форма функції розподілу, але для 27 випадків, коли важливі декілька механізмів розсіювання, для обчислення цих величин необхідно використовувати масове моделювання Монте-Карло. Можна виразити потік енергії, який з’являється в рівнянні через тензор температури. Підводячи підсумок, рівняння балансу для функції дрейфованого максвеллівського розподілу спрощуються до 1 n n n S t q      ( 1.4) 0 1 (( ) ) ( )e d n z m W T W nkT k E W W t x x               ( 1.5) 2 22 1 ( )n d i n m q nq nm nkT E t m x m              ( 1.6) У цих виразах 20,5 1,5dW nm nkT  n de n  1.3.3 Метод Монте – Карло Результати розрахунків коротких уніполярних пристроїв показують, що легування витоку/відіграє важливу роль з точки зору струму приладу, що є насамперед ефектом послідовного опору. З представлених результатів очевидно, що нестаціонарний транспорт відіграє меншу роль у FD SOI довжиною 90 нм пристроїв, тоді як важливість нестаціонарного транспорту та пов’язане з ним перевищення швидкості різко зростає для пристрою FD SOI із довжиною затвора 14 нм. Ці результати свідчать про те, що необхідно включити рівняння енергетичного балансу, якщо потрібно досягти належного моделювання нанорозмірних пристроїв із довжиною затвора менше 100 нм. 28 Ще одним питанням, яке заслуговує на подальшу увагу, є залежність результатів моделювання від вибору часу релаксації енергії. Ми бачимо сильну залежність поточного струму від вибору часу релаксації енергії для найменшої структури, що моделюється, що свідчить про необхідність правильного визначення часу релаксації енергії. Час релаксації енергії, у свою чергу, є параметром, що залежить від зміщення та геометрії, і його точне визначення неможливо. Неможливість правильно визначити час релаксації енергії в моделях гідродинамічного/енергетичного балансу було основною мотивацією для розробки моделей на основі частинок. У попередньому розділі ми розглянули континуальні методи опису транспорту в напівпровідниках, зокрема моделі дрейфу-дифузії та гідродинамічні моделі, які виводяться з моментів напівкласичного рівняння переносу Больцмана (BTE). Як наближення до BTE, було показано, що у випадку малих пристроїв такі підходи стають неточними або повністю виходять з ладу. Дійсно, можна уявити, що, оскільки фізичні розміри зменшуються, на певному рівні опис континууму струму порушується, і дискретна природа окремих частинок заряду, що становить щільність заряду в області активного пристрою, стає важливою. Мікроскопічне моделювання руху окремих частинок у присутності сил, що діють на них через зовнішні поля, а також внутрішні поля кристалічної решітки та інші заряди в системі вже давно популярне в науковому співтоваристві, де молекулярна динаміка моделювання атомів і молекул вже давно використовується для дослідження термодинамічних властивостей рідин і газів. У твердих тілах, таких як напівпровідники і Відомо, що в металах транспорт переважає випадкове розсіювання через домішки, коливання гратки тощо, які рандомізують імпульс і енергію заряджених частинок у часі. Отже, стохастичні методи моделювання цих випадкових подій розсіювання особливо корисні для опису транспорту в напівпровідниках, зокрема метод Монте-Карло. Багаточасинковий метод Монте-Карло використовуються досить довго, як чисельний метод для моделювання нерівноважного транспорту в напівпровідникових матеріалах і 29 пристроях, а також був предметом численних книг і оглядів [36]. У застосуванні до транспортних задач випадкове блукання генерується за допомогою алгоритмів генерації випадкових чисел, звичайних для сучасних комп’ютерів, для моделювання стохастичного руху частинок, що піддаються процеси зіткнення. Цей процес генерації випадкових блукань є частиною дуже загальної методики, яка використовується для оцінки інтегральних рівнянь, і пов’язана із загальною технікою випадкової вибірки, яка використовується для оцінки багатовимірних інтегралів [37]. Основна техніка, застосована до проблем транспортування, полягає в моделюванні вільного руху частинок (відомого як вільний політ), що завершується миттєвими випадковими подіями розсіювання. Алгоритм Монте-Карло полягає в генеруванні випадкового часу вільного польоту для кожної частинки, вибір типу розсіювання, що відбувається в кінці вільного польоту, зміна кінцевої енергії та імпульсу частинки після кетерінг, а потім повторення процедури для наступного вільного польоту. Вибірка руху частинок у різні моменти часу протягом моделювання дозволяє статистичну оцінку фізично цікавих величин, таких як функція розподілу окремої частинки, середня швидкість дрейфу в присутності прикладеного електричного поля, середня енергія частинок тощо. імітація ансамблю частинок, Отже, для дослідження фізичної ситуації, яка викликає інтерсе, може бути змодельована нестаціонарна залежна від часу еволюція розподілу електронів і дірок під впливом рушійної сили, що залежить від часу. 30 2. МОДЕЛЮВАННЯ РОБОТИ ВАРІЗОННОГО ДІОДА НА ОСНОВІ InPAs 2.1 Модель діода 2.1.1 Опис зонної структури та параметрів напівровідника Ключовою характеристикою матеріалу є його закон дисперсії. Від виначає як саме пов’язана кінетична енергія електрона з хвиотрвим вектором. У випадку, коли енергія електрона далоек авід енергії за якої виникає пробій у напівпровіднику, цю залежність можна описати досить модельно. Такий підхід зображено на рис. 2.1 Рис. 2.1. Аналітична модель для опису залежності кіненичної енергії носіів заряду від хвильового числа, яка вкирисанга при моделюванні Зона провідності за такого наближення описується трьохдолинною ізотропною 6 6 6 C C CL X   – моделлю зони провідності, Наьбільщ пийнаятна аналітична залежність ( )k може бути представлена як на Так як багаточастинковий метод Монте-Карло використовує [44]: 2 2 ( )(1 ( )) 2 k k k m   , (2.1) 31 де  – є так званий коефіцієнт непараболічності. В найжчій долині зони провідності його можна описати у вигляді 2(1 ) g e m E m     , де gE – є величина ширина забороненої зони. Скадні потрійні сполуки 1z zA B C , насправді поєднують властивості двох бінарних сполук AC та BC , які відповідають значенням 0z  та 1z  . Виходячи з цих міркувань можна визначити більшість партері потірйних сполук. Основним законом, який встановлює залежність від складу параметрів матеріалу є правилом Вегарта( правило пропорційного розподілу). Тоді, коефіцієнт непараболічності 1z zA B C  :   1 1 z zA B C AC BCz z       , (2.2) енергія мінімумів долин:     1, , ,1 1 z zg A B C g AC g BC вигE zE z E z z C       , (2.3) де вигC – емпіричний параметр, що визначається виходячи з експериментуж; постійна решітки:   1 (1 ) z zA B C AC BCa x za z a     , (2.4) В подібний спосіб вихначеється і ефектинна маса елетронів * * * 1 1 ( ) AC BC z z m z m m    (2.5) Для дірок закон дисперсії описуюється анізотропнми законом дисперсії [46]: 2 2 2 2 4 2 2 2 2 2 2 0 1 2 3 2 0 ( ( ) 4 12( )( )) 2 V y zk k k k k k k k m           E E , (2.6) де 1 , 2 и 3 – є параметрами Латтінжера. Іх визначення для сполуки AzB1- zC в залежності від z відьбуваться через парамктри 1 , 2 и 3 [43], які визначаються з правила Верагта. 32 Більщ складну залежність від від z використовується для опису густини сполуки     1 1 3 3 3 1 z z z z AC AC BC BC A B C A B C z a z a x a         , (2.7) та статичної і високочастотної діелектричних сталих: 1 1 1 1 2 1 z z z z z z A B C A B C A B C F F        , 1 (1 ) z zA B C AC BCF zF z F     , (2.8) де: 1 2 AC AC AC F      ; 1 2 BC BC BC F      . 2.1.2 Вразуовання розсіяння носіїв заряду Зазвичай вважається, що зовнішні сили, що діють на носії заряду пов’язані з кристалічною решіткою, малі, а, отже, можна застосувати теорію збурень [24], згідно з якою ймовірність переходу зі стану з хвильовим вектором i k в стан з хвольовим вектором f k , спричинене дією j - го механизму визначається золотим правилом Фермі: 2 22 2 ( , ) ( ) ( ) i if i f f i j f ik k W k k M f V iE E E E , (2.9) де i ifk k M - матричний элемент переходу, j V - потенциал збуорення, який обумовлений дією дефекту, i E - энергія системи до розсіяння, f E - энергія системи після розсіяння. Для фононної система, яка перебуває у стані термодинамічної рівноваги число фононів визначається планковим розподілом: 1(exp( / ) 1) qb qb N kT (2.10) В такому разі, імовіоність переходу зі стану в стан за рахунок розсіювання на фононах: 33 2 22 1 1 ( , ) ( ) 2 2 2 qb if i f f iqb qb C W I k k N NM E E (2.11) Отже, для визначення імовірності переходу электрона з початкового стану k , якому він має енергію k E у будь-який інший стан k з енергією k E . який відбувається з поглинанням або випроміненням фонона з енергією qb треба визначити параметр звязку qb C для даного механізму та інтеграл перекриття 2( , )I k k , а в кінці обчислити інтеграл от(2.11) по усім можливим кінцевим станам k : 2 22 1 1 ( ) ( , ) ( ) 2 2 2 qb i qb k k qb qb C W k I k k N dk NM E E , (2.12) знак вказує на процесс поглинаня фонона, а знак на процесс випромінювання. В роботі розглянуто: полярне оптичне росіяння, розміяння на акустичних фононах, розсіяння на деформаціонному потенціалі оптичних фононів в L – долині та внутридолиннета міджолинне росзіяння; розсіяння на домішках; розсіяння на сплавному потенціалі. Відоповідно до (2.12) Іміовіність розсіяння деформаціному потенціалі акустичних фононів буде мати вигляд: 2 2 2 2 2 2 1 1 1 ( ) ( ) 4 2 2 2 d ДА qb qb l m q W k N q dq s k m                          , (2.13) Тут q – модуль хвильового вектору фонона; d – ефективна повздовжня константа деформаційного потенціалу; ls – швидкість повздовжніх акустичних хвиль;  – густина. Якщо акустичні фонони збуджуються в потрійній сполуці 1z zA B C , то необхідну величину швидкості акустичної хвилі можна знайти за формулою[45]: 34      1 1 22 1 1 1z z z z AC BC A B C A B C AC BC z A H z z D z C H a zH z H           , (2.14) Коефіцієнтами, що входять до виразу визначаються за формулами: , 2 2 , , вч BC BC LO BC ст BC H     ; , 2 2 , , вч AC AC LO AC ст AC H     (2.15) BC AC C D B AH H B     , де AC AC A a        , 3 AC AC BC BC a B a             , 2 BC BC C a        . (2.16) Межі інтегрування в (2.13) визначають кут, на який розсіюється частика: 22 (1 2 ) 1 cos 1 2 S lq q m s k          , (2.17) де: * /S lq m s . Для розсіяння на деформаційному потенціалі оптичних фононів ймовірність розсіяння набуває вигляду:     0 2 *3 2 0 3 1 1 1 2 22 d ДО D m W N K k                 , (2.18) де: 0D – величина деформаційного потенціалу;  – частота оптичних фононів; * dm – ефективна маса густини станів в зоні, в яку відбувається розсіяння; Коефіцієнт 1 0,5 1 ( ) 1 2 1 2 K                   враховує непараболічність зони та дорівнює одиниці у випадку коли вона параболічна. Значення кута, на який розсіюється електрон знаходитися в межах:  0 02 1 2 / 1 cos 1 2 im q q k             . (2.19) Якщо оптичні фонони збуджуються в потрійній сполуці 1z zA B C , то велчина енергії, яку буду мати оптичні фонони визначається за формулою[45]: 35  1 1 1 1 2 , , ,2 2 2 2 , , , , , , 1z z z z z z ст A B C вч AC вч BC LO A B C LO AC LO BC вч A B C ст AC ст BC z z                          , (2.20) де ст , вч –відповідно статична та високочастотна діелектричні сталі. У слабких електичних полях в полярних напівпровідниках домінуючим механізмом є розсіяння на коливаннях протилежно заряджених атомів кристалу. Таке розсіяння отримало назву полярне оптичне розсіяння. Імовірність того, що електрон буде розсіяним у цьому випадку є :     2 * 2 2 max 2 * min 1 1 1 1 ln 2 2 ПО e m q kk W k N k q m                                 ст ,(2.21) Параметри k  визначаються співвідношеннями     2 1dk m        , а кут, на який буже розсіяний електрон, знаходиться за формулою (2.19) Окремо розглядається взаємодія з короткохвильовими фононами на межі зони Брилюена, для яких частота не залежить від хвилевого вектору. Ймовірність такого розсіяння обяислюється аналологічно до (2.18) у вигляді:        3 2 2 2 1 11 1 1 1 2 2 1 22 i i j jj ij di МД i j j j j j Z D m W N                   , (2.22) У виразі врахована кількість еквівалентних долин у зоні Брілюена jZ в які можливий перехід електронів. Всі переходи відбуваються із виконная закона збереження енергії j i ij    , де велчина ij є енергетична відстань між i та j , між яким відбувається перехід. У потрійних сполуках AzB1-zC виникає невпорядкованість розміщення атомів А та В у вузлах кристалічної решітки, що є причною порушення періодичності кристалічного потенціалу. З цим пов’язана поява додаткового механізму розсіяння - розсіяння на сплаві [49]. Імовірність розсіяння електрона визнається на підставі моделі моделі Мотта: 36       3 2 0 * * 4 1 ( ) 1 2 2 1СП СП x x r U d W m m               , (2.23) де СПU – потенціал розсіювання на сплаві, 0r – радіус комірки Вігнера- Зейтца. Емпіричний параметр d називається невпорядкованістю кристалічної решітки і вибирається з міркувань співпадіння результатів розрахунку та експериментальних данний для потрійних сполук, зокрема (0 1d  ). В умовах високої концентрації носіів заряду, а відповідно і високим ступенем легування актульаним механізмом розсіяння є розсіяння на заряджених домішках. В роботі для визначення імовірності розсіяння електрона використана формула Рідлі [43]:         2 1 exp 2 m ВН ДР m k b W k W k b k               , (2.24) де 1/3 max 0,5 npb N  ; dN – концентрація домішок, а  ВНW k – імовірність розсіяння, яка визначається за наближенням Брукса-Херрінга: . 2.1.3 Алгоритм методу Монте-Карло На початку моделювання визначаються початкові значення енергії та імпульсів частинок. Імпульс частинки визначається як: 1 2 3 2 ( 1,5) T p p z z z , (2.25) де 1 2 3 , ,z z z - випадкові числа, що відповідають рівномірному розподілу в інтервалі від 0 дло 1, щ дає змогу отримати величину, яка з точністю до 10 %, має нормальний розподіл, T p m kT Частинка має кінетичну енергією  t і рухається під дією зовнішніх сил до моменту розсіяння, який буде визначатися імовірністю розсіяння 37   W t , де   W t є сумою ймовірностей усіх механізмів розсіяння. Момент часу t , в який відбувається розсіяння визначається умовою: 1 0 ( ( )) ) ln t W t dt zE (2.26) Пока нерівність (2.26) не виконується частинка рухається у просторі відповідно до рівнянь:    * 2 * 2 1 ee e e e p tdr dt m p t m    , (2.27)     * * * 1 2 e e ee e e e m E mdp E e x dt m E              , (2.28)    * * ( ) p p G p p dp E e x E x m dt m         . (2.29)  p p p dr E p dt   , (2.30) Рівняння (2.27) – (2.30) враховують координатну залежність параметрів матеріала. Градієнт складу ( )z x визначає градієнт усіх інших величин:  x , ( )GE x ,  CE x ,  VE x , ( )x . Розподіл електростатичного потенціалу знаходиться як розвязок рівняння Пуассона:   ( , )x y      (2.31) Для цього застосовуєтсья повний багатосітковий метод[ 40]. Енергія та імпульс частинки, які вона має в момент розсіяння визначають імовірність сцьогд процесу і вибір конкретного механізму розсіяння. Для цього генерується випадкового число z , яке рівномірно розподілене у інтервалі [0,1] та його значення порівнюється з відповідними сумами: 38 1 1 1 ( ) ( ) ( ) ( ) m m i i i i W W z W W        , (2.32) де max1,2,3...m m , maxm – кількість механізмів розсіювання. При виконні умови (2.31) обирається механізм з номером m . Після цього відповідно до вибраного механізму знаходиться нові значення імпульсу та енергії, які частикна набуде в результаті розсіяння. 2.2. Результати моделювання 2.2.1 Постановка зададачі Для визначення максимальних частот коливань досліджено діоди довжиною від 0,5 до 1,2 мкм. Діоди являють собою структуру типу n+-n- n+-. За даними [13], концентрація в активній n-області становила 10 17 см -3 . Області n+- є сильно легованими областями поблизу контактів з концентрацією донорів 2⸱10 17 см -3 . Важливість анодної концентрації була продемонстрована в [13]. Тому в основній частині розглянутих діодів концентрація 2⸱10 17 см -3 була використана для досягнення високої величини сили, що діє на отвори в напрямку анода. Координатна залежність складу сплаву (мольна частка миш’яку) z в InP1- zAsz представлена таким чином:     0 2 0 02 0, 1 exp( ), 2 x x z x x x x x           , (2.33) де 0x - координата гетероінтерфейсу між однорідним InP і варізонному InPAs. Показано, що оцінка ефективності генерації є максимальною для складу As рівною 0,8 на контакті анода. З цих причин ми формуємо складений розподіл, щоб отримати значення z= 0.8 на аноді шлязом вибору  . Розподіл (2.33) дещо схожий на гетероперехід, особливо у випадку діода 39 довжиною 500 нм. З цієї причини очікується, що збільшення енергії електрона в активній області та перенесення в долини супутника будуть досить швидкими). Моделювання діодів здійснюється за допомогою методики ансамблю Монте-Карло (EMC). Основна перевага цього методу полягає в його здатності враховувати високочастотні процеси в діоді та швидкі зміни прикладеної зовнішньої напруги в терагерцовому діапазоні. Також для цих цілей використовується метод синхронного ансамблю для визначення часу вільного прольоту носіїв заряду. Зона провідності представлена трьома долинами непараболічності (Г, L, X). Залежність кінетичної енергії електрона від модуля хвильового вектора описується з використання (2.1)   2 2 1 ( ) 2 ( ) K K k E z E m z     , (2.34) де ( )m z - ефективна маса електрона, ( )z коефіцієнт непараболічності та , – модифікована постійна Планка. Тут ( )m z та ( )z є функціями складу матеріaлу. Залежності потенціальної енергії від складу були взяті у формі:   21.35 1.094 0.1ГE z z z   (2.35)   2.05 0.97LE z z  (2.36)   2.21 0.84XE z z  (2.37) Оскільки мінімальна енергія, яка повинен мати електрон для здійснення акту ударної іонізації thE визначається дірками з найбільшою ефективнобю масою, то валентна зони представлена саме підзоною важких дірок. Зонна діаграма діода довжиною 500 нм наведена на рис.1. 40 Рис. 2.2. – Зоннна діаграма діода, на прикладі діода довжиною 0,5 мкм. Відомо, що можливість використання напівпровідника з вузькою забороненою зоною в діоді Ганна обмежує ударну іонізацію в області сильного поля. У цьому випадку густина струму буде зростати з часом і через деякий час осциляції Ганна можуть стати некогерентними [39]. В InAs електрони можуть нагріватися електричним полем до енергій, що значно перевищують енергію забороненої зони. Навпаки, в InAs n-типу швидкість нагрівання дірок обмежена через квазіелектричне поле та відносно пласку зону важких дірок. Як свідчить рис. 2,2 використання однієї зони для опису валентної зони можливе, якщо на діод буде прикладена напруга зміщення. Порогова енергія ударної іонізації в InPAs залежить від кристалографічної орієнтації та обмежена фактором непараболічності [40].       21 ( ) 1 1 4 ( ) ( ) 2 / 1 . 2 ( ) th gE z z E z z                , (2.38) де ( )z  це відношення ефективної маси електрона до ефективної маси дірки. Імовірність іонізації ( )IIW E визначається приблизно таким чином 0( , ) ( )( ( ))n II thW E z W z E E z  , (2.39) 41 де 0( )W z і n залежать від матеріалу та відомі як параметри прилягання або м’якості. Підхід Келдиша ( 2n  )використовується, що вимагає жорсткого порогу для ударної іонізації. У цьому звіті розглядається початковий етап ударної іоназації, і це наближення можна застосувати. Кристалічна орієнтація (напрямок від катода до анода) пропонується [1,1,1], в результаті ударна іонізація розглядається лише в Г-долині та Х-долині. Розглянуто такі найважливіші механізми розсіювання, як розсіювання на деформаційному потенціалі неполярних оптичних фононів та акустичних фононів, розсіювання полярних оптичних фононів, міждолинне розсіювання, розсіювання сплаву та іонізованих домішків. Параметри матеріалу, використані в чисельній моделі, відповідають [12]. Для симуляції діода була використана 2D модель. Покращена схема розподілу заряду «хмара в комірці» (CIC) була реалізована для підвищення просторової точності за наявності просторово залежної діелектричної проникності разом із зменшенням власної сили. Для визначення 2D- розподілу потенціалу використовувався повний багатосітковий (FMG) метод розв’язування рівняння Пуассона, що забезпечує максимально швидке отримання результатів [40]. 2.2.2 Характеристки діодів на постійному струмі Залежність щільності струму від прикладеного зміщення для діодів на основі InP/InPAs без (w/o) II та з II, для різних довжин діодів 500 нм, 640 нм та 1280 нм, показано на рис. 2. 42 Рис.2.3. – Current density J versus bias voltage U for diodes with difference length, with II – 1-3, w/o II – 4-6;. 1,4 – 500 nm, 2,5 – 640 nm, 3,6 – 1280 nm. Зі спостережуваного збільшення щільності струму можна зробити висновок, що II є домінуючим фактором для формування I-V – характеристики. У розглянутому діоді ударна іонізація ініціюється тільки одним типом носія. Це характерно відрізняється від випадку, коли обидва типи носіїв можуть викликати ударну іонізацію. Це призводить до експоненціального зростання струму без лавинного пробою. Енергетична діаграма гетероструктурного діода InP/InPAs, енергетичний розподіл носіїв зміни та кількість подій ударної іонізації (в умовних одиницях) при напрузі зміщення 1 В, наведено на рис.2.4. 43 Рис.2.4 – Енергетична діаграма діода, розподіли носіїв заряду по долинах і кількість ударних актів іонізації в діоді довжиною 500 нм при напрузі зміщення 1 В. Залежності ( )CE z та ( )VE z визначати основні компоненти квазіелектричних полів у зоні провідності та валентній зоні відповідно. Як видно на рис. 3, іонізація відбулася поблизу контакту анода. У цій області квазіелектричне поле у валентній зоні діє в протилежному напрямку порівняно з іншими частинами діода. Електрони та дірки тягнуться в одному напрямку до анода. Менша частина дірки рухається до катода, але не може отримати достатньо кінетичної енергії для здійснення акту ударної іонізації. Таким чином, дірковий струм є алгебраїчною сумою струму, пов’язаного з дірками, що рухаються в анод, і струму через дірки, які пройшли до катода. Оскільки на процес ударної іонізації впливає енергія носія, ймовірність ударної іонізації зменшується зі збільшенням довжини діода. Це пояснюється як збільшенням частоти розсіювання фононів, так і меншою силою, що діє на електрони. В результаті наростання струму за рахунок ударної іонізації в діоді 1280 нм не спостерігається. Причому в цьому випадку при великій напрузі зміщення відбувається невелике падіння струму. У коротших діодах диференціальний опір постійного струму додатний, і струм зростає зі 44 зменшенням довжини діода. Лінійне зростання струму в цьому масштабі зумовлене домінуванням ударної іонізації лише одного типу носіїв (електронів). Насправді коефіцієнт розмноження Mn для електронів не може перевищувати 2, і ймовірність того, що дірка зробить хоча б один акт ударної іонізації, практично дорівнює нулю. Таким чином, залежності, отримані без врахуванням ударної іонізації можна розглядати як спосіб енергетичної релаксації. 2.2.3 Частотні характеристики діодів Енергетичні та частотні характеристики діодів можна оцінити, розглянувши їх роботу в режимі коливань. У цьому режимі передбачається, що діодна структура розміщена в резонаторі, а прикладена напруга U(t), що діє на структуру, є сумою постійної напруги зміщення U0 та альтернативної напруги U(t). Пропонується, щоб резонатор мав високу добротність Q, тому U(t) являє собою першу гармоніку з частотою f і амплітудою U1: 0 1( ) sin2U t U U ft (2.40) Через відносно сильний струм канавок при високій напрузі зміщення діапазон напруги, що відповідає режиму тривалості, може бути обмежений зверху. В основному це має місце для короткого діода, який завжди стабільний при U1= 0. Нестабільність струму зростає при зміщенні, яке перевищує деяке критичне значення. оцінюється понад 1,25 В для діода 500 нм, і понад 1,3 В для діода 500 нм, і це залежить від частоти. Розподіл електричного поля вздовж діода 500 нм у різничні моменти часу протягом періоду коливань показано на рис.2.5. 45 Рис.2.5. – Розподіл напруженості електричного поля в діоді з довжиною 640 нм, U0= 2.25 V, U1=1 V: 1 – t =0; 2 – t = T/4; 3 – t = T/2; 4 – t = 3Т/4. Розподіли отримано при оптимальних значеннях напруги та частоти резонатора. Основним видом нестабільності струму в діодах є утворення шарів накопичення. Шари накопичення утворюються в основному в другій половині градієнтного шару, де напруженість електричного поля вища, особливо в анодній області 200 нм. Ця ширина відповідає області, де домінуючим процесом є ударна іонізація (II). Проте в діоді 1280 нм автоколивання струму спостерігаються і при напрузі вище 2,6 В. Цей факт можна пов'язати зі статичним негативним опором у цих діодах. Частота коливань близько 150 ГГц. При меншому зміщенні коливання напруги затухають. У режимі резонансу для діода 1280 нм перевищує 1,3 В. Високочастотні властивості діода оцінювали шляхом визначення ефективності коливань ( ) на відповідній частоті резонатора. розраховувався як відношення потужності, що генерується на частоті резонатора (перша гармоніка), до середньої потужності постійного струму. Щоб знайти максимальну ефективність, була проведена оптимізація шляхом зміни зміщення U0 та амплітуди U1. Для оцінки впливу ударної іонізації на 46 генерацію діода також отримано характеристику без (без) ударної іонізації. Залежності оптимізованої ефективності коливань від частоти резонатора наведено на рис. 5 та рис. 6. Концентрація донорів у всіх діодах однакова і дорівнює 10 17 см -3 . Рис.2.6. - Частотна залежність максимальної ефективності коливань (1,2) та густини потужності змінного струму (1’,2’) в діоді 1280 нм: 1,1’- залежності, отримані із врахуванням ударної іонізації I, 2,2’- залежності, отримані без врахуванням ударної іонізації. Як видно з рисунку, піковий ККД для довгого діода може досягати 10 % і відповідає f0=100 ГГц. Коливання спостерігаються в діапазоні від 70 до 250 ГГц. Оптимальна напруга зсуву для всіх розглянутих діодів становить (2...3) Up, що відповідає ефекту Ганна в однорідному діоді. Частота коливань зростає зі зменшенням довжини діода, f0 =180 ГГц для діода 640 нм і f0 = 250 ГГц для діода 500 нм. Максимальний ККД 3 % і 1,3 % відповідно. Діоди мають широкий діапазон частот. Так, для діода 640 ГГц вона становить від 120 до 350 ГГц, для діода 500 нм – від 140 до 400 ГГц. Найвище значення граничної частоти суми розглянутої довжини можна отримати для діода 500 нм і становить 400 ГГц. Це значення вище, 47 ніж у градієнтних діодах на основі InGaAs з такими ж параметрами, розглянутими в [13], для яких гранична частота становила 350 ГГц. Рис.2.7. Частотна залежність максимального ККД коливань (1-4) і густини потужності змінного струму (1’-4’) в діоді довжиною 500 нм (1,1’,2,2’) та 640 нм (3,3’,4,4’) відповідно: 1,1’,3,3’ –залежності, отримані із врахуванням ударної іонізації I, 2,2’,4,4’- залежності, отримані без врахуванням ударної іонізації Очікується, що II впливатиме на верхню межу частоти коливань. Насправді підвищення граничної частоти в порівнянні з діодами без ударної іонізації спостерігається тільки в довгих діодах. Таким чином, частота зрізу була досягнута до 250 ГГц (200 ГГц у випадку імітаційних діодів без ударної іонізації). У коротких діодах гранична частота залишається такою ж, але величина ККД нижче, ніж при моделюванні діодів без ІІ. Як видно з результатів моделювання, наведених на рис.2.6 і рис.2.7, ККД коливань в діодах з ударної іонізації зсувається в бік менших значень у всіх розглянутих випадках. У той же час напруга зміщення, що відповідає максимальному ККД в діоді з ударною іонізацією нижче, ніж у діода без ударної іонізації. 48 Розглянуто вплив ударної іонізації на можливість отримання генерації. У деяких випадках цей процес може призвести до поліпшення частотних властивостей діода. Наприклад, це має місце в градуйованому діоді на основі InGaAs [13]. У нашому розгляді використовується проста модель II, яка враховує процес релаксації лише в нижній Г – долині та сателітній Х-долині. Фактично існує різка анізотропія ударної іонізації. Тому ударна іонізація в L – долині також може відігравати важливу роль у релаксації енергії гарячих електронів. Цей випадок необхідно проаналізувати в майбутньому. Це може призвести до зміни наших оцінок. Другим важливим фактором, що впливає на результати, може бути ефект саморозігріву, що виникає внаслідок сильного поля в анодній області. Передбачається, що це буде важлива роль у довгому діоді, і його також слід враховувати при моделюванні. 49 ВИСНОВКИ Представлено результат моделювання методом Монте-Карло діодів InPAs з варізонним проміжком Наше дослідження показало, що варизонні діоди з InP на катоді та вузькозонний напівпровідник InP0.2As0.8 на анодному контакті можуть бути використані як ефективні активні елементи для генерації коливань у нижній частині терагерцового діапазону. Ефектом переносного електрона вважається резон виникнення коливань. Гранична частота генерації архівів 400 ГГц, як діод довжиною 500 нм. Показано, що діоди мають широкий діапазон частот понад 100 ГГц. Розглянуто вплив ударної іонізації на можливість отримання генерації. У деяких випадках цей процес може призвести до поліпшення частотних властивостей діода. Наприклад, це має місце в градуйованому діоді на основі InGaAs ПЕРЕЛІК ПОСИЛАНЬ 50 1. L-F. Shi, A. Zahid, A. Ren, M.Z. Ali, H. Yue, M.A. Imran, Y. Shi, & Abbasi. The perspectives and trends of THz technology in material research for future communication - a comprehensive review // Physica Scripta. 2023. Vol. 98, No. 6, P. 065006-1 - 065006-22. 2. M. Asada, S. Suzuk. Terahertz Emitter Using Resonant-Tunneling Diode and Applications // Sensors. 2021, Vol. 21, No. 4, P. 1384-1 – 1384-19. 3. R. Izumi, S. Suzuki, M Asada. 1.98 THz resonant-tunneling-diode oscillator with reduced conduction loss by thick antenna electrode // 42nd International Conference on Infrared, Millimeter, and Terahertz Waves (IRMMW-THz 2017), art. no. 17259158, 1 (Cancun: IEEE: 2017). 4. M.S. Vitiello, A. Tredicucci. Physics and technology of Terahertz quantum cascade lasers // Advances in Physics: X. 2021. Vol. 6, No. 1. P. 1893809-1 -1893809- 31. 5. A. Khalatpour, A.K. Paulsen, C. Deimert, Z.R. Wasilewski, Q. Hu. High- power portable terahertz laser systems // Nat. Photonics. 2021. Vol. 15, P. 16-20. 6. A. S. Hajo, O. Yilmazoglu, A. Dadgar, F. Küppers and T. Kusserow. Reliable GaN-Based THz Gunn Diodes With Side-Contact and Field-Plate Technologies // IEEE Access. 2020. Vol. 8, P. 84116-84122. 7. Capasso, F. (1988). Compositionally Graded Semiconductors and their Device Applications. In: Margaritondo, G. (eds) Electronic Structure of Semiconductor Heterojunctions. Perspectives in Condensed Matter Physics, vol 1. Springer, Dordrecht. 8. I. Storozhenko. Gunn Diodes Based on Graded-Gap GaInPAs // Journal of Nano- and Electronic Physics. 2020. Vol. 12, No1. P. 01015-1 – 01015-9. 9. Prykhodko K. H., Zozulia V. О., Botsula O. V. Graded band gap ingaas diodes for terahertz applications // YSF–2017: 2017 IEEE International Young Scientists Forum on Applied Physics and Engineering, October 17- 20, 2017: Conference Paper. Lviv, Ukraine, 2017. P. 291–294. 51 10. O. V. Botsula, K. H. Prykhodko. Graded Band InGaN- Based Diode for Noise Generation in Terahertz Range // 2020 IEEE Ukrainian Microwave Week, UkrMW 2020, September 21–25, 2020: Conference Paper. Kharkiv, Ukraine, 2020. P. 925–928. 11. I. Storozhenko, S. Sanin. Advanced Micron Sized Gunn Diode Based on Graded-Gap GaPAs // GaInAs. Journal of Nano- and Electronic Physics. 2022. Vol. 14, no. 1. P. 01027–1–01027–5. 12. O.V. Botsula, K.H. Prykhodko, V.A. Zozulia. InGaAs- based graded gap active elements with static cathode domain for terahertz range // Journal of Nano- and Electronic Physics. 2019. Vol. 11, No1. P. 01006-1–01006-5. 13. O. Botsula, K. Prykhodko, V. Zozulia. Impact Ionization in Graded Gap Transferred Electron Diode // 2021 IEEE 3rd Ukraine Conference on Electrical and Computer Engineering, (UKRCON 2021). August 26–28, 2021: Conference Paper. Lviv, Ukraine, 2021. P. 504–508. 14. S. Asmontas, V. Gružinskis, P. Shiktorov, et al. Monte Carlo simulation of THz frequency range Gunn effect in InP MOSFET at impact ionization conditions // The 31st International Conference on the Physics of Semiconductors (ICPS 2012), 2012. P. 385 – 386. 15. G. M. Dunn, G. J. Rees and J. P. R. David Monte Carlo simulation of impact ionization and current multiplication in short GaAs - diodes// Semiconductor Science and Technology. 1997. Vol. 12, No. 1. P. 11 – 120. 16. G. M. Dunn, A. Philips and P. J. Topham Gunn instabilities in power HEMTs// Semicond. Sci Technol.2001. Vol. 16, No. 1. P. 262 – 566. 17. H. Marinchio et al., Experimental and theoretical investigation of terahertz optical-beating detection by plasma waves in high electron mobility transistors Phys. Stat. Sol. (c). 2008 Vol. 5, No. 1. P. 257 – 259. 18. K. Brennan, K. Hess, High field transport in GaAs, InP and InAs// Solid- State Electron.1984. Vol. 27, No. 4. P. 347 –354. 19. A.Dargys, J. Kundrotas, Handbook on Physical Propertie of Ge, Si, GaAs and InP, Vilnius, Sience and Encyklopedia Publishers (1994). 52 20. V. Gružinskis, E. Starikov, P. Shiktorov. Monte Carlo simulation of THz frequency Gunn effect in GaAs MOSFET with excess of electrons in channel at impact ionization conditions // Lithuanian Journal of Physics. 2014. Vol. 54, No. 1. P. 7 – 10. 21. J.H. Zhao et al, Monte Carlo Simulation of 4H-SiC IMPATT Diodes Semicond. Sci. Technol. 2000. Vol. 15, P. 1093–1100. 22. F. Osaka, T. Mikawa and T. Kaneda, "Impact ionization coefficients of electrons and holes in ,// IEEE Journal of Quantum Electronics. 1985, Vol. 21, No. 9. P. 1326-1338. 23. D.K. Ferry, S.M. Goodnick, J. Bird. Transport in Nanostructures, // 2009. – Cambridge University Press; 2nd edition, – 2009, 670 P. 24. D. Vasileska, S.M. Goodnick. Computational electronics // Materials Science and Engineering, Reports: A Review: Journal. 2002, R38. No. 5, P. 181 – 236. 25. D. Vasileska and S. M. Goodnick. Computational Electronics // Morgan and Claypool, - 2006, 216 P. 26. A. Schuetz, S. Selberherr, H.W. Poetzl. A Two-Dimensional Model of the Avalanche Effect in MOS Transistors // Solid-State Electron. 1982, Vol. 25. P.177-183. 27. P. Antognetti and G. Massobrio. Semiconductor Device Modeling with SPICE // McGraw-Hill, New York, - 1993, 479 P. 28. M. Shur, Physics of Semiconductor Devices (Prentice Hall Series in Solid State Physical Electronics) // Prentice-Hall, - 1990, 704 P. 29. D. L. Scharfetter and D. L. Gummel. Large-signal analysis of a silicon Read diode oscillator // IEEE Transaction on Electron Devices. 1969, Vol. 16. No. 1, P. 64 - 74. 30. R. Lake, G. Klimeck, R.C. Bowen, and D. Jovanovic. Single and multiband modeling of quantum electron transport through layered semiconductor devices // Journal of Applied Physics. 1997, Vol. 81. No. 12, P. 7845 – 7869. https://www.amazon.com/s/ref=dp_byline_sr_book_3?ie=UTF8&field-author=Jonathan+Bird&text=Jonathan+Bird&sort=relevancerank&search-alias=books 53 31. G. Baccarani and M.R. Wordeman. An Investigation of Steady-State Velocity Overshoot in Silicon // Solid-State Electronics. 1985, Vol. 28. No. 4, P. 407 – 416. 32. S. Cordier. Hyperbolicity of grad’s extension of hydrodynamic models for ionospheric plasma part i: the single species case // Mathematical Models and Methods in Applied Sciences. 1994, Vol. 4. No. 5. P. 625 - 645. 33. K. Tomizawa. Numerical Simulation of Submicron Semiconductor Devices // The Artech House Materials Science Library, - 1993, 341 P. 34. H. K. Gummel. A self-consistent iterative scheme for one-dimensional steady state transistor calculations // IEEE Transactions on Electron Devices. 1964, Vol. 11. No. 10. P. 455 - 465. 35. T. M. Apostol. Calculus, Vol. II, Multi-Variable Calculus and Linear Algebra // Blaisdell, Waltham, - 1969, 673 P. 36. T. Grasser, T.-W. Tang, H. Kosina, S. Selberherr. A review of hydrodynamic and energy-transport models for semiconductor device simulation // Proceedings of the IEEE. 2003, Vol. 91. No. 2. P. 251 – 274. 37. C. Jacoboni and P. Lugli. The Monte Carlo Method for Semiconductor Device Simulation // Springer-Verlag, Vienna, - 1989, 359 P 38. Vurgaftman I., Meyer J. R., Ram-Mohan L. R. Band parameters for III–V compound semiconductors and their alloys // Jornal of applied physics. 2001. Vol. 89, No11. P. 5815 – 5875. 39. Levinshtein M., Kostamovaara J., Vainshtein S. Breakdown phenomena in semiconductoes and semiconductors devises Selected topics of electronic and system // Singapore: World Scientific Publishing, 2005. 208 p. 40. Brennan K., Mansour N. Monte Carlo calculation of electron impact ionization in bulk InAs and HgCdTe // Jornal of applied physics. 1991. Vol. 69, No. 11. P. 7844–7847. 41. Joppich W., Mijalkovic S. Multigrid Methods for Process simulation // New York: Springer-Verlag, 1993. 313 p. 54 42. Pérez, S., González, T., Pardo, D., Matéos, J. (2008). Terahertz Gunn-like oscillations in InGaAs/InAlAs planar diodes. Journal of Applied Physics, 103(9). 43. Jacoboni, C., Reggiani, L. (1983). The Monte Carlo method for the solution of charge transport in semiconductors with applications to covalent materials. Reviews of Modern Physics, 55(3), 645–705. 44. Lundstrom, M. (2000). Fundamentals of Carrier Transport. 45. Nederveen, K. (1989). Ensemble Monte Carlo simulation of electron transport in AlGaAs/GaAs heterostructures. PhDT. 46. Vurgaftman, I., Meyer, J. R., Ram-Mohan, L. R. (2001). Band parameters for III–V compound semiconductors and their alloys. Journal of Applied Physics, 89(11), 5815–5875 47. Hockney, R., Eastwood, J. (1988). Computer simulation using particles. CRC Press. 48. Littlejohn, M. A., Hauser, J. R., Glisson, T. H., Ferry, D. K., Harrison, J. W. (1978). Alloy scattering and high field transport in ternary and quaternary III–V semiconductors. Solid-state Electronics, 21(1), 107–114. 49. Harrison, J. W., Hauser, J. R. (1976). Alloy scattering in ternary III-V compounds. Physical Review, 13(12), 5347–5350.