УДК 519.684, 004.42

МОДЕЛИРОВАНИЕ КИНЕТИЧЕСКИХ СВОЙСТВ ДВУМЕРНЫХ ПОЛУПРОВОДНИКОВЫХ МАТЕРИАЛОВ с ПРОИЗВОЛЬНОЙ ФОРМОЙ ПЕРВОЙ ЗОНЫ БРИЛЛЮЭНА МЕТОДОМ МОНТЕ-КАРЛО

, ,

Волгоградский государственный технический университет

*****@***ru

Разработана программа, позволяющая на основе квазиклассического моделирования методом Монте-Карло вычислять плотность постоянного тока и среднее время релаксации в двумерных материалах в условиях воздействия постоянных и переменных электрических и магнитных полей, использующая для распараллеливания расчетов сопроцессор Intel Xeon Phi.

Ключевые слова: метод Монте-Карло, кинетические свойства материалов, Intel Xeon Phi.

Исследование кинетических свойств низкоразмерных материалов позволяет предсказывать новые эффекты и определять перспективы их применения в электронных устройствах. Наряду с расчетами на основе решения кинетического уравнения для этого используется имитационное моделирование методом Монте-Карло [1-3], которое в квазиклассической ситуации состоит в следующем.

1) Задаем начальное распределение квазиимпульсов набора частиц исходя из равновесной функции распределения [1].

2) На каждом шаге по времени численно решаем квазиклассические уравнения движения частицы и вычисляем вероятность её рассеяния с новым значением квазиимпульса.

3) Зная вероятность рассеяния частицы, разыгрываем момент столкновения носителя заряда с нерегулярностями решетки и, исходя из законов сохранения энергии и квазиимпульса, вычисляем новые значения компонент квазиимпульса и скорость.

НЕ нашли? Не то? Что вы ищете?

4) После прохождения достаточно большого числа шагов по времени усредняем значение скорости по всему временному промежутку.

5) Повторив расчет для большого числа частиц, усредняем скорость по ансамблю.

6) Зная среднее значение скорости, вычисляем постоянную составляющую плотности электрического тока.

Особенностью метода является возможность проведения расчетов для каждой из частиц в отдельном потоке, что позволяет распараллелить задачу на большое число процессорных ядер. Ситуация осложняется тем, что часто вероятность рассеяния задается довольно сложной формулой, а вызов функции для её расчета происходит сравнительно часто. Выходом может явиться предварительный расчет вероятности рассеяния и последующая интерполяция по уже известному массиву вероятностей. При таком подходе наиболее эффективно [4] использование сопроцессоров Xeon Phi, обладающих большим числом процессорных ядер и общей памятью, что позволяет распараллелить расчеты при минимальном изменении последовательного алгоритма.

Рисунок 1 ‑ Геометрия задачи

Рассмотрим вычисление вольтамперной характеристики (ВАХ) полупроводниковой сверхрешетки (CР). Будем рассматривать систему координат , оси которой повернуты относительно кристаллографических осей структуры на угол против часовой стрелки (рисунок 1, слева). Первая зона Бриллюэна в этом случае представляется прямоугольником ABCD (рисунок 1, справа), а энергетический спектр структуры имеет вид:

(1)

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

Вычисление вероятностей рассеяния носителей заряда на фононах сводится к нахождению интеграла следующего вида:

. (2)

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

. (3)

Используя свойства дельта-функции, преобразуем (3):

, (4)

где ‑ i-й корень уравнения

(5)

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

а)

б)

Рисунок 2 ‑ Сечение графика энергии плоскостью равной энергии

Проверка правильности расчетов вероятности проводилась для случая в (1), для которого известно аналитическое выражение для полной вероятности рассеяния. В случае акустических фононов

. (6)

Здесь - константа деформационного потенциала акустических фононов, - поверхностная плотность материала, - скорость звука, ‑ абсолютная температура, . В случае оптических фононов

(7)

Здесь ‑ константа деформационного потенциала оптических фононов, ‑ энергия оптических фононов. Интеграл, входящий в (6), (7), берется точно:

(8)

Здесь - эллиптический интеграл первого рода, - полный эллиптический интеграл первого рода. Трехмерные графики зависимости вероятности рассеяния электронов на акустических и оптических фононах от компонент квазиимпульса приведены на рисунках 3 а) и 3 б), соответственно ().

а) на акустических фононах

б) на оптических фононах

Рисунок 3 ‑ Полная вероятность рассеяния электрона

Зависимость вероятности рассеяния на оптических фононах имеет характерную область , которая соответствует значениям энергии электрона, меньшим энергии фонона. На обоих графиках вероятность имеет разрывы второго рода. Математически это связано с равенством нулю знаменателя в (4) при определенных значениях импульса. При достижении вероятностью значения, соответствующего достоверному событию (в случае нормированной вероятности это значение равно 1) при расчете считаем, что происходит рассеяние. Вычисление средней скорости (по времени и ансамблю) проводилось по методике, аналогичной использованной в [2, 3].

При расчете ВАХ СР в (1) полагали , вектор напряженности электрического поля направлен вдоль оси СР, г, эВ, эВ/см, эВ, см, эВ, , г/см2, см/с. Время измеряем в единицах c, скорость измеряем в единицах см/с. В этом случае константы (в единицах ), компоненты скорости

, (9)

. Длина промежутка времени, в течение которого следим за движением электрона, составляет 2000, шаг по времени 0.0005. В результате моделирования получается зависимость плотности постоянного тока вдоль оси СР от напряженности постоянного продольного электрического поля, представленная на рис. 4 а). Компонента плотности тока измеряется здесь в некоторых условных единицах, напряженность электрического поля ‑ в единицах (здесь ед. СГСЭ – заряд электрона). На рисунке 4 б) представлена зависимость среднего времени свободного пробега от напряженности электрического поля. Известно (см., например, [5]), что положение максимума вольтамперной характеристики подчиняется условию

(10)

(в приближении постоянного времени релаксации ). В нашем случае (приблизительно 6 ед. СГСЭ), и левая часть соотношения (10) составляет , а правая, полагая в соответствии с графиком на рисунке 4 б) , равна , то есть (10) довольно точно выполняется.

а) ВАХ СР

б) Среднее время релаксации

Рисунок 5 ‑ Результаты моделирования

Таким образом, реализована программа, которая позволяет численно находить полную вероятность рассеяния носителей заряда на фононах в случае материала со сложной формой первой зоны Бриллюэна, проведено сравнение результатов работы этой программы с расчетами по явной формуле. Во-вторых, реализован метод Монте-Карло для расчета плотности постоянного электрического тока с использованием распараллеливания расчетов на сопроцессоре Intel Xeon Phi при помощи технологии OpenMP, проведена проверка программы на примере вычисления вольтамперной характеристики полупроводниковой свехрешетки. Проведенные проверки указывают на соответствие результатов основным физическим закономерностям, программа правильно работает для известных предельных случаев и может быть применена для моделирования кинетических явлений в сложных конфигурациях электрических и магнитных полей.

Работа выполнена при финансовой поддержке Минобрнауки России на выполнение государственных работ в сфере научной деятельности в рамках базовой части государственного задания № 000/411 (коды проектов: 522 и 3154) с использованием сопроцессора Intel Xeon Phi, приобретенного в рамках Программы стратегического развития Волгоградского государственного технического университета.

Литература

1.  Кашурников, В. А., Красавин, методы квантовой статистики. – М.: ФИЗМАТЛИТ, 2010. – 628 c.

2.  Завьялов, моделирование эффекта выпрямления тока, индуцированного электромагнитной волной в графене / , , // ФТП. – 2010. – Т. 44. – Вып. 4. – С. 910-914.

3.  Konchenkov, V. I. Influence of constant electric field on circular photogalvanic effect in material with Rashba Hamiltonian / V. I. Konchenkov, S. V. Kryuchkov, D. V. Zav’yalov // J Comput Electron. – 2014. – V. 13. – P. 996 – 1009.

4.  Rezaur Rahman, Intel Xeon Phi Coprocessor Architecture and Tools – NY.: Apress Open, 2013 – 220c.

5.  Басс, свойства полупроводников со сверхрешетками / , , . ‑ М. : Наука, 1989. – 288 с.

MONTE CARLO MODELING OF KINETIC PROPERTIES OF TWO-DIMENSIONAL SEMICONDUCTOR MATERIALS WITH ARBITRARY SHAPE OF Brillouin zone

R. I. Gushchin, D. V. Zav`yalov, V. I. Konchenkov, M. V. Skachkov

Volgograd State Technical University

*****@***ru

The program is developed which can calculate on the base of Monte Carlo method a current density and a mean relaxation time in two-dimensional materials under the influence of constant and alternative electric and magnetic fields using paralleling on coprocessor Intel Xeon Phi

Keywords: Monte Carlo simulations, kinetic properties of materials, Intel Xeon Phi.