Горячие электроны в нитриде индия: новый метод численного решения задачи электронного транспорта Текст научной статьи по специальности «Физика»

ФИЗИКА КОНДЕНСИРОВАННОГО СОСТОЯНИЯ ВЕЩЕСТВА

Горячие электроны в нитриде индия: новый метод численного решения задачи электронного транспорта

H.A. Масюков0, A.B. Дмитриев0

Московский государственный университет имени М. В. Ломоносова, физический факультет, кафедра физики низких температур и сверхпроводимости. Россия, 119991, Москва, Ленинские горы, д. 1, стр. 2.

E-mail: 11 masyukov@mig.phys. msu. ru, b dmitriev@lt.phys. msu.su

Статья поступила 28.11.2008, подписана в печать 10.04.2009.

Проведен численный расчет электронного транспорта в InN в греющем электрическом поле. Полученные результаты очень хорошо согласуются с известными экспериментальными данными. Для вычислений использован новый численный метод решения транспортного уравнения Больцмана в изотропном пространственно-однородном случае, основанный на аппроксимации функции распределения из точек относительно редкой сетки в импульсном пространстве и итерационном поиске стационарного решения.

Ключевые слова: кинетика, электронный транспорт, горячие электроны, нитрид индия, численные расчеты.

УДК: 538.935. PACS: 72.20.Ht.

Введение

Исследование с решеткой типа вюрцита стало особенно актуально после того, как в работах [1] было предложено новое значение ширины запрещенной зоны ег ~ 0.7 эВ, а не 2 эВ, как считалось ранее. Это открывает перспективы широкого применения сплавов нитрида индия и широкозонного нитрида галлия в оптоэлектронике. Соответственно возрастает интерес к исследованию электронных свойств и, в частности, к изучению транспортных свойств этого полупроводника в сильных электрических полях. Эта проблема уже привлекала внимание исследователей: так, в работах [2-7] численно решалась задача нелинейного электронного транспорта в методом Монте-Карло для относительно небольших концентраций электронов N = 1017 см^3. В то же время имеющиеся экспериментальные данные были получены в работе [8] для плотностей носителей, почти на два порядка больших. Соответственно ни в каких расчетных работах нет сравнения полученных численно результатов с экспериментальными данными. К тому же часть этих работ [2, 3] основана на старых представлениях об как о широкозонном по-

лупроводнике. Изменение такой фундаментальной величины, как ширина щели, и уточнение других параметров спектра, осуществленное в последнее время, могут оказать существенное влияние на вид расчетных вольт-амперных характеристик и существенно приблизить их к экспериментальным результатам, как и использование адекватной величины концентрации электронов. Именно такие расчеты и проведены в настоящей работе.

В ней предложен и применен к модели массивного образца 1п]Ч и-типа вычислительно эффективный итерационный метод решения уравнения Больцмана для пространственно-однородного изотропного случая. Полученные численно результаты сравниваются с экспериментальными данными [8].

1. Постановка задачи и физическая модель

Наборы параметров 1пГ\Г, использованные в работах [2-7, 9], представлены в табл. 1. В наших расчетах были учтены все эти наборы параметров и проведено сравнение полученных результатов.

В соответствии с экспериментальной ситуацией [8] концентрация носителей заряда была взята равной

Таблица 1

Параметры материала

Величина Символ Единица измерения Значение

HI [9] Н2 [2, 3] НЗ [4, 6, 7]

Плотность Р г/см3 6.81

Скорость звука S 105 см/с 2.53 3.78

Статическая диэлектрическая проницаемость Хо 6.7 8.4 8.4

Высокочастотная диэлектрическая проницаемость Хоо 9.3 15.3 15.3

Эффективная масса т* т0 0.07 0.12 0.045

Запрещенная зона % эВ 0.7 2.0 0.8

Деформационный потенциал Н эВ 3.6 7.1 7.1

Пьезомодуль ец К/и2 0.375

Энергия оптического фонона hwо эВ 0.073

Расстояние до ближайшего побочного минимума зоны проводимости А эВ > 1.77

N = 9- 1018 см^3. Отметим, что при такой концентрации электронный газ остается вырожденным во всем рассматриваемом интервале температур меньше 300 К. Концентрация заряженных примесей считается равной концентрации электронов.

Поскольку, по имеющимся представлениям, побочные экстремумы зоны проводимости расположены высоко (см. табл. 1), задача решается в однодолинном приближении (справедливость его более подробно обсуждается в разделе 3). Закон дисперсии электронов взят изотропный и непараболический:

пЧ2

2т*

; е( \ + ае),

п дк + \т)соП

с условием нормировки функции распределения ]Г/(Й) = АП/.

(1)

(2)

(

экранирования, ц=±(Ы —к), и¡(ц) — закон дисперсии фононов: оо(ц) = зц для БА и РА, оо(ц) = щ для РО; Щ = Ыч(оо(ц), Т) — число фононов с волновым вектором ц при температуре решетки Т. Разогрев решетки не принимается во внимание и считается, что температура решетки остается постоянной и не зависит от электронной температуры Те. Для фононных механизмов рассеяния верхний и нижний знаки относятся к испусканию и поглощению фонона соответственно. Обычно рассеяние на акустических фононах можно считать упругим. В этом случае БА- и РА-рассеяния могут быть рассмотрены как один механизм с вероятностью рассеяния

¥-ар (к^к')

где е — энергия электрона, й — его волновой вектор, фактор а учитывает непараболичность спектра и зависит от ширины запрещенной зоны и эффективной массы т* на дне долины: а = (1 — т*/пц)2

Для описания электронного транспорта мы используем кинетическое уравнение Больцмана. В пространственно-однородном случае это интегро-дифференциаль-ное уравнение имеет вид

1

V а~

(В0А(?)+Вра(?))(2Лг, + 1)х

6(е(к)-е№).

Здесь Е — приложенное к образцу электрическое поле, /(¿, й) — функция распределения электронов, в — заряд электрона, V — объем образца,

%) „ = "►*)/(*')(!-/(*))-

А' -»*')/(*)(! — /(Л'»] (3)

есть интеграл столкновений, Ж/(к к') — вероятность перехода электрона из состояния к в к'.

В наших расчетах учитываются все основные механизмы рассеяния электронов в 1пГ\Г: рассеяние на заряженных примесях (г), на деформационных акустических (БА), поляризационных акустических (РА) и поляризационных оптических (РО) фононах. В табл. 2 представлены вероятности рассеяния для рассматриваемых механизмов [10] в том виде, в котором они будут использованы далее. Мы пользовались следующими обозначениями: х= (Хж~Хо ') '> = 47ге14/хо, А — длина

Интеграл столкновений (3) не линеен по функции распределения, что существенно в случае вырожденной статистики. В этом случае методы, основанные на сведении уравнения Больцмана к системе линейных уравнений на основе линейной по функции распределения аппроксимации интеграла столкновений, становятся неприменимыми. Метод Монте-Карло, хоть и может быть обобщен на случай вырожденной статистики, делается вычислительно громоздким. Стоит также отметить, что линеаризация интеграла столкновений с рассматриваемыми вероятностями рассеяния в принципе не может быть проведена без дополнительных предположений, так как вероятности рассеяния зависят от длины экранирования А, которая в свою очередь зависит от функции распределения:

-1/2

хит

4тге2

Хо

Е

д№

(4)

Таким образом, известные методы не вполне удобны для решения поставленной задачи. Предлагаемый в работе метод требует лишь гладкости функции распределения.

2. Интеграл столкновений

Прежде чем перейти к рассмотрению использованного численного метода, необходимо преобразовать интеграл столкновений (3), учитывая аксиальную симметрию функции распределения ¡(к)=[(к^,к±). Запишем интеграл столкновений в виде

дt)m

д[\ 0'Лр

Вероятности рассеяния

Таблица 2

Механизм рассеяния \М(к^к')\2 В(Я)

г ЫУ2^\М(к^к") 2х х 6(е(к) - е(к')) ( 1 4тге2 ц2 \2 1 и хо Я2 + \-2) ф

БА х 6(е(к) - е(к') т М?)) Ж^ + А-О ЙН2

РА Нер)2 , 2 Р8 4

РО Ншоа^2 У

и рассмотрим каждый член отдельно в сферической системе координат.

Для первого слагаемого имеем

2тг

2тт ( Aire2

й V Хо

где введены вспомогательные функции

ЫО(кЩ(Щ,к±)11ф(Щ,к±),

m*k

dk'k'2(E(k') -e(k))

4ir3h

2añ2k2

m"

Ц(Ь\\,Ь±)--

Ü(khk±):

2тг

2lT

d eos e(f(k eos в, k sin в) -f(k||, k±)), d<p

(2к2^2Щк cos 6^2k±k sin в cos ф+Х^2) 2k2 - 2k\\k cos в + Л^2

3/2 '

({2к2^2Щк cos 9+\^2)2^(2k±k sin в)2)

причем k = + k2±. Таким образом, в члене (df/dt)¡

остается интегрирование только по одной переменной в конечных пределах. Далее,

/<9/\ _ 2тт / /<9/\ absorP rdf \dt)ро ñ \\dt/po \dt

of y

PO

(

л? absorp,

d/Vmis = D(k'±)I*0±(khk±)I*0±(khk±),

dtJ po I¡°+(k hk±)

d eos 9[f(k'+eos 9, k'+ sin0)(l ^f(klbk±))(N0+ 1) - /(A,,, k±)(1 - /(k'+ eos в, k'+ sin 0))Л*>],

/PO±

d eos в [/ (k'_ eos в, k'_ sin в) (1 - / (А,,, k±))N0 --f(k||,A±)(1 —f{k'_ cosd,k'_ sin0))(iVo+ 1)],

2

2тг

(k¡¡tk±):

d4>Bp0(q±)

q±

q±2 + X-2

где ц± = д±(к,к'±,в,ф) = (к2 + к - 2к\1к'±совв -— 2к±к'± Б'твсовф)^2, к'± определяются из уравнений е(к±) = е(к) ± йшо, Л^о — число оптических фононов в кристалле. Нетрудно заметить, что в члене (9//с?0ро остается интегрирование по двум переменным в конечных пределах.

Для акустических фононов в упругом приближении

где

1вкР(к п.

d cos 0(/ (к cos в, к sin в) - f (к» ,k1_)),

lf(khk±)-

dф{Bш{q)+Bш{q)){2Nq + \)

Kq2 + X^2J ' о

q = q(k, в, ф) = (2k2 - 2кф cos в - 2k±k sin в cos ф)1/2.

Таким образом, и в члене (df/dt)^Р остается интегрирование по двум переменным в конечных пределах. Заметим, что в том случае, когда для рассеяния на акустических фононах недостаточно ограничиться упругим приближением, в соответствующем члене становится необходимо выполнять интегрирование по трем переменным. Это существенно повышает вычислительную сложность метода. Однако неупругость DA-рассеяния может быть легко учтена приближенно, если пренебречь ее влиянием на величину импульса фонона, что обусловлено малостью его скорости по сравнению со скоростью электрона, следующим образом. Используя тогда ту же зависимость q = q(k, в, ф), что и в упругом приближении, мы интегрируем по двум изоэнергетическим поверхностям с

8F

к'±(к,в,ф)ск±(-

ñsq(k, в, ф).

Таким образом, в полной аналогии с РО-рассеянием, рассмотренным выше, в члене (9//с?0ар удается получить интегрирование только по двум переменным в конечных пределах.

3. Численный метод

Задача нахождения функции распределения решается на двумерной сетке в импульсном про-

странстве к'„

+ (i + 1/2) Д0)

/До,

г = 0,1,..., 2ктах/Д0 - 1, / = 0,1,... ,ктах/А0, а размер ктах области, где задана сетка, выбирается так, чтобы функция распределения полностью помещалась внутри этой области. Разумеется, далекие экспоненциально малые «хвосты» распределения неизбежно окажутся за пределами области размером ктах, но число описываемых ими носителей и их роль в транспорте тоже экспоненциально мала, так что ими можно пренебречь. Исключение могут составлять задачи типа ионизации, где важны именно электроны с большой энергией, но здесь мы такими проблемами не занимаемся. Значения функции распределения между точками сетки находятся с помощью бикубической аппроксимации. Отметим, что при гладкой функции распределения и хорошей ее аппроксимации между точками сетки сама эта сетка может быть довольно редкой.

Интересующее нас стационарное решение уравнения Больцмана находится с помощью итерационного метода, построенного на основе зависящего от времени уравнения (1). На итерации можно с этой точки зрения смотреть как на временную эволюцию функции распределения в соответствии с уравнением (1).

В момент времени t = 0 задается некоторое начальное приближение = ,к'±). Итерационный шаг с п-го на (и+1)-й временной слой организован следующим образом. Функция распределения интерполируется из узлов сетки в (к^,к±)-пространство, для этого используется бикубическая сплайн-интерполяция. Интерполированные значения функции распределения

¡п(крк±) подставляются в формулу (4) для вычисления Л" = Х(1п(к^,к±)) и в интеграл столкновений. Значения

интеграла столкновений вычисляются только

для узлов сетки по формулам раздела 2. Заметим, что интегрирование может быть проведено с необходимым шагом, никак не связанным с шагом сетки. Переход на следующий временной слой осуществляется вычислением новых значений функции распределения в узлах сетки непосредственно из транспортного уравнения Больцма-на (1)

Vj, 107см с

! 1 _

еЕ

И~

At, кj

-At

coll

и перенормировкой (2). Итерационная процедура повторяется до достижения сходимости.

В качестве начального приближения мы для простоты использовали фермиевскую функцию Е(к,Т,[1 о) с решеточной температурой и химическим потенциалом, соответствующим нулевой температуре, умноженную на коэффициент, определяемый из условия нормировки (2).

Для рассмотренного диапазона электрических полей Е мы использовали ктах ~ 1.2 • 107 см^1. Функция распределения при этом хорошо локализована в области симуляции (рис. 1). Однодолинное приближение оправдано тем, что е(ктах) <С А для всех наборов параметров материала, для которых мы проводим моделирование.

к 107см-1

1. Функция распределения. " с, isin, = 5-Ю-13 с,

параметры материала HI

0.6 -0.8 к±, 10 см"

Сетка 15 х 30,

Рис.

At = Ю-15 с! isim = 5-"КГ1"3 с, Т = 77 К, Е = 15 кВ/ем,

Сходимость метода проверялась в следующих тестах (рис. 2). При фиксированных значениях поля и температуры решетки Е = 15 kB/cm, Т = 77 Кис набором параметров материала HI (табл. 1) проводились вычисления с различными At и числом узлов сетки. Оптимальными оказались шаг At = 1(Г15 с и сетка размером 15 х 30. В этом случае время симуляции 4im = Ю^12 с соответствует icpt_: — 2 мин компьютерного времени (Intel Core 2 Duo Т72500, 2 ГГц, использовалось одно ядро). Взяв в качестве начального приближения неравновесную функцию, например Еп/Щ(к) = F(k,T,^о) + cos(10^6ft), можно было убедиться, что результат не зависит от начального приближения (рис. 2). На рис. 2 также можно видеть, что результат слабо зависит и от того, в упругом

- -J—___4____ _ 4- - - -j >

— 15x30, At = le- 15 с: tCPU= 107 с □ 20x40, At = 1 е - 15с: tCPU = 226 с + 15x30, At = 2е - 15 с: tCPU= 53 с

О 15x30, At = 0.5е - 15 с: tCPU = 218 с

- - 15x30, At = le - 15 с, q/el: tCPU= 511 с .........15x30, At = le - 15 с, n/eq: tCPU= 111 с

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

0.9 1.0

xlO"12

Рис. 2. Иллюстрация сходимости. Зависимость дрейфовой скорости иа от времени симуляции 4|т для различных параметров численного метода. Параметры материала Н1, Т = 77 К, Е = 15 кВ/см

или неупругом (ц/е1) приближении было учтено рассеяние на акустических фононах (в последнем случае имеется в виду приближенный учет неупругости, описанный в разделе 2). Таким образом, акустическое рассеяние может быть описано в рамках упругого приближения, и мы так и делали в дальнейших расчетах.

4. Некоторые результаты

Для нахождения полевой зависимости дрейфовой скорости были использованы различные наборы параметров материала (табл. 1). Из сравнения результатов

vd, 107см с 1

Те> К

Е, кВ см

Рис. 3. Полевая зависимость дрейфовой скорости (а) и электронной температуры (б). Результаты вычислений для различных наборов параметров материала и экспериментальные данные. Сетка 15x30, = 10-ь с, 4т1 = 5 • Ю-13 с, Г = 77 К

вычислений с экспериментальными данными [8] видно (рис. 3,о), что набор Н2 дает наилучшее соответствие. Стоит заметить, что набор Н2 включает старые (до того как работы [1] были опубликованы) значения параметров зонной структуры, которые были получены на относительно «грязных» образцах с концентрациями электронов, близкими к используемой в наших вычислениях величине А^ = 9 • 18 см^3. На рис. 3,6 представлено сравнение расчетной полевой зависимость электронной температуры Те при температуре решетки Т = 77 К с результатами, полученными в работе [8] на основе анализа экспериментальных данных. Способ нахождения Те описан в Приложении.

На рис. 4 представлена вычисленная полевая зависимость дрейфовой скорости при различных температурах решетки для набора параметров Н1. Был использован именно этот набор параметров, потому что он соответствует современным представлениям о зонной структуре ¡пГ* [11].

уа, 107см с-1

хЮ

+ 4К □ 150 К ^300 К

Й о

2 4 6 8 10 12 14 16 Е, кВ см-1

Рис. 4. Половая зависимость дрейфовой скорости для различных температур решетки. Сетка 15 х 30, Д* = Ю-15 с, 4|ш = 5 • Ю-13 с, Н1

Заключение

Проведен численный расчет электронного транспорта в в греющем электрическом поле. Полученные результаты очень хорошо согласуются с известными экспериментальными данными. Для вычислений использован новый численный метод решения транспортного уравнения Больцмана в изотропном пространственно-однородном случае, основанный на аппроксимации функции распределения из узлов относительно редкой сетки в импульсном пространстве и итерационном поиске стационарного решения.

В работе метод применен к с концентрацией электронов N = 9- 1018 см^3. Определены оптимальные параметры численного метода (размер сетки, шаг итераций и т.п.). Рассчитаны полевые зависимости дрейфовой скорости для различных температур решетки. Были проведены вычисления с наиболее часто используемыми наборами параметров материала. Наилучшее соответствие результатов расчетов с экспериментальными данными [8] получается при использовании набора Н2 (табл. 1).

Приложение

Электронная температура Те и химический потенциал ц. определяются путем аппроксимации истинной функции распределения смещенной фермиев-

ской функцией Р(к,Те,//.). Это делается следующим образом. Пусть мы определили среднее значение волнового вектора (йц) на п-м шаге итерационной процедуры. Тогда можно определить среднюю энергию теплового движения электронов V. Рассмотрим функции

ё\(Теф)

1

7T2N

dkk2F(k,Te,ß), (5)

dkk2e(k)F(k,Te,n). (6)

Требуется найти такие значения Те и ц., что

gi(Te,ß) = l, g2(Te,ß) = U.

Выберем в качестве начального приближения некоторые Те0 и //о, такие, что gi(Te0, ц0) = Q, g2(Teo,/io) = i/o-Раскладывая (5) и (6) no Те — TeQ и [1 — ßQ, получим систему линейных уравнений для нахождения Те и ц.

(Sgl Sgl ;>т ф. Ofc Öfö ОТ, ;>,,

Производные берутся в точке (Тео,//о), каждая из них содержит одно интегрирование. Эта система легко решается, решение сходится за несколько итераций при хорошем начальном приближении. В качестве начального приближения естественно выбрать значения искомых величин, найденные на предыдущем шаге главной итерационной процедуры.

Список литературы

1. Dauydou V.Yu., Klochikhin A.A., Seisyan R.P. et al. // Phys. Stat Sol. (b). 2002. 229, N 1; 2002. 230, N 4.

2. Bellotti E., Doshi B.K., Brennan K.F. et al. // J. Appl. Phys. 1999. 85, N 2. P. 916.

3. Foutz B.E., O'Leary S.K., Shur Al.S., Eastman L.F. // J. Appl. Phys! 1999. 85, N 11. P. 7727.

4. O'Leary S.K., Foutz B.E., Shur Al.S., Eastman L.F. // Appl. Phys. Lett. 2005. 87. P. 222103.

5. Starikov E., Shiktorov P., Gruinskis V. et al. // J. Appl. Phys. 2005. 98. P. 08370L

6. Polyakov V.M., Schwierz F. // Appl. Phys. Lett. 2006. 88. P. 032101.

7. Polyakov V.M., Schmerz F. // J. Appl. Phys. 2006. 99. P. 113705.

8. Zanato D., Balkan N.. Ridley B.K. et al. // Semicond. Sei. Technol. 2004. 19. P. 1024.

9. Hsu L., Jones P.E., Li S.X. et al. // J. Appl. Phys. 2007. 102. P. 073705.

10. Гашпмахер В.Ф., Левинсон И.Б. Рассеяние носителей тока в металлах и полупроводниках. М., 1984.

11. Walukiewicz W., Auger J.W., Yu. K.M. et al. // J. Phys. D. 2006. 39. P. 83.

Hot electrons in wurtzite indium nitride: a novel numerical method for solving the problem of electron transport

N.A. Masyukov", A.V. Dmitriev

Department of Low-Temperature Physics and Superconductivity, Faculty of Physics, M. V. Lomonosov Moscow State University, Moscow 119991, Russia.

E-mail: 11 masyukov@mig.phys. msu. ru, b dmitriev@lt.phys. msu.su.

Hot electron transport in InN is studied using a novel numerical procedure based on the interpolation of the electron distribution function from the grid points in the momentum space and on the iterative solution of the Boltzman transport equation. The agreement with the existing experimental data is excellent.

Keywords: kinetics, electron transport, hot electrons, indium nitride, numerical calculations.

PACS: 72.20.Ht.

Received 28 November 2008.

English version: Moscow University Physics Bulletin 4(2009).

Сведения об авторах

1. Масюков Никита Андреевич — студент; e-mail: masyukov@tnig.phys.msu.ru.

2. Дмитриев Алексей Владимирович — докт. физ.-мат. наук, профессор; тел. 939-59-05, e-mail: dmitriev@lt.phys.msu.su.