МОДЕЛИРОВАНИЕ ПОВЕДЕНИЯ ГРАФЕНА ВО ВНЕШНИХ ЭЛЕКТРИЧЕСКИХ ПОЛЯХ

Панферов Анатолий Дмитриевич; Новиков Николай Андреевич; Трунов Александр Алексеевич

ББК 32.972+22.379 ГРНТИ 27.35.57 УДК 519.688:004.942

А. Д. Панферов, Н. А. Новиков, А. А. Трунов

Моделирование поведения графена во внешних электрических полях

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

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

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

Ключевые слова и фразы: численное моделирование, графен, функция распределения

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

Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 18-07-00778.

Введение

Интерес к графену, как одному из перспективных материалов для современной электроники, в значительной степени обусловлен особенностями его отклика на действие внешних электрических полей [1-4]. Для его практического использования необходимо уметь достаточно точно соотносить параметры действующего поля (частотный спектр, напряженность, зависимость от времени и т.п.) с результатами его действия на этот материал. Это требует адекватной теоретической модели и методов её использования. На основе классической модели безмассовых фермионов [5,6] в работах [7-9] был предложен и развит подход, основанный на использовании квантового кинетического уравнения для описания на языке квантовой теории поля процессов образования и эволюции носителей заряда в псевдо-двумерном монослое графена в условиях действия внешнего классического электрического поля. В результате для функции распределения носителей f (р1,р2,1) была определена в явном виде система обыкновенных дифференциальных уравнений (см., например, [7]) с выражающимися через параметры действующего поля переменными коэффициентами. Эта функция распределения определяется в двумерном импульсном пространстве —п < р1 < п, —п < р2 < п (первая зона Бриллюэна) и характеризует вероятность обнаружить носитель заряда со значениями импульса {р1, р2} в момент времени Зная эту функцию, можно вычислять наблюдаемые параметры, такие как поверхностные плотности носителей заряда и тока.

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

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

В предлагаемой работе показывается, что характер зависимости I(Р1, Р2) во время действия внешнего электрического поля существенно отличается от I(Р1,Р2,£ ^ то). В присутствии внешнего поля значения I(Р1,Р2) доступны для вычисления с использованием стандартной арифметики двойной точности во всей области определения —п < р1 < п, —п < р2 < п или в достаточно большой её части. При этом I(Р1,Р2) ведет себя гладко за исключением области максимальных значений. Такое поведение функции распределения позволяет определять местоположение области максимальных значений для использования в ней плотных регулярных сеток или адаптивной сетки. Отдельно исследован вопрос об особенностях моделирования действия слабых электрических полей. Представлены результаты реализации процедуры воспроизведения всей эволюции рассматриваемой модели из начального состояния I(Р1,Р2,£ ^ —то) = 0 в конечное состояние I(р1, р2, ^ ^ то) = 0 и определения по этим данным динамики наблюдаемых характеристик.

1. Функция распределения носителей в импульсном пространстве в присутствии внешнего поля

Исходно в качестве цели использования кинетической модели [7] графена рассматривалось воспроизведение его состояния после завершения действия внешнего поля по аналогии с приложением этого подхода к квантовой электродинамике [11,12]. Но при сопоставлении полученных результатов с данными экспериментальных исследований и результатами использования других теоретических моделей стало понятно, что протекающие во время действия поля процессы также крайне интересны и существенны [13-15]. Это стало поводом рассматривать время £ как полноценный третий аргумент функции распределения и исследовать её поведение на всем интервале эволюции. С точки зрения проблемы выбора расчетных сеток в импульсном пространстве существенным оказалось то, что в присутствии даже относительно слабого внешнего поля характер поведения функции распределения претерпевает существенные изменения.

Рассмотрим в качестве примера моделирование результатов действия короткого импульса электрического поля вида

(1)

/ \

Е^) = Ео сов(ш*)ехр( - у^у), Е2(*)=0,

где циклическая частота ш = определяется через частоту действующего поля V, параметр а определяет характерную длительность импульса, а амплитудное значение Ео соответствует максимальной напряженности поля в момент времени £ = 0. Поле имеет постоянное направление действия, с которым ассоциировано направление первой координатной оси.

Выберем следующие значения параметров: длительность ~ 2 х 10-12 с, частота несущей ~ 2 х 1012 Гц, амплитуда ~ 3 х 104 В/см. Это достаточно характерный набор параметров, определяющий

1(Р1,Р2)

1

0.01 0.0001 1x10' 1x10' 1x10"

1x10' 1x10'

,-14

Рисунок 1. Функция распределения f (р1,р2,Ь ^ — то) после завершения действия импульса поля с параметрами: длительность ~ 2 х 10-12 с, частота несущей ~ 2 х 1012 Гц, амплитуда ~ 3 х 104 В/см

очень короткий высокочастотный импульс с напряжённостью поля в окрестностях границы нелинейности. Результаты действия на графен

импульсов с близкими параметрами экспериментально исследовались, например, в работе [2]. По завершению такого импульса функция распределения принимает вид I(Р1,Р2,£ ^ то), представленный на рисунке 1. Для наглядности на этом рисунке на регулярной сетке размером 312 х 312 показана вся области определения I(р1,р2). Не смотря на то, что логарифмическая шкала имеет диапазон значений 1.0 х 10-14 < I < 1.0, только в 62 точках (из представленных на рисунке 97344) значение функции распределения превышает порог отсечения 1.0 х 10-14.

Более детальное представление о характере поведения I (р1,р2) в центральной области импульсного пространства даёт рисунок 2. На нем на сетке с размерами 102 х 102 показана область локализации значений I(р1,р2), определяющих наблюдаемые характеристики модели.

1(Р1,Р2)

Рисунок2. Детализация функции распределения f (р1,р2,£ ^ -<х), представленной на рисунке 1

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

области определения мы не просто не знаем значений этой функции, но и не имеем возможности их вычислить при условии ограничения использованием переменных типа ¿впЫв на 64-х разрядных вычислительных системах. Можно выйти за рамки этих ограничений и использовать "длинную арифметику" с произвольной точностью представления операндов [16]. Но это ведёт к неприемлемому для практического использования снижению быстродействия [17]. Поэтому исходно был реализован поиск областей локализации с последовательным перебором регулярных сеток увеличивающейся плотности [10].

Теперь рассмотрим поведение функции распределения в присутствии электрического поля. Для определенности и наглядности выберем момент времени £ = 0, когда напряженность электрического поля (1) максимальна. При выбранных ранее параметрах I(р1,р2,£ = 0) имеет вид, представленный на рисунке 3.

1(Р1,Р2)

1 г

0.01 Г

Рисунок 3. Функция распределения f (р1,р2,1 = 0) для момента времени с максимальным значением напряжённости внешнего поля. Остальные параметры аналогичны рисунку 1

Для полноты картины на рисунке 4 представлен аналогичный "предшественник" для рисунка 2.

1(Р1,Р2)

Рисунок4. Детализация функции распределения f (р1,р2,4 = 0), представленной на рисунке 3, для центральной области импульсного пространства

Хотя область локализации относительно больших значений, близких или сопоставимых с максимально достижимым /(Р1,Р2) ^ 1.0, остаётся очень компактной, вне неё значения функции не уходят за пределы нижней границы 1.0 х 10-14 доступного для вычислений с использованием переменных типа ¿впЫв диапазона. При этом / (^1,^2) ведет себя гладко.

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

2. Моделирование действия слабых внешних полей

Представленные результаты демонстрируют возможность эффективно локализовать и точно воспроизводить функцию распределения носителей, формируемую импульсом поля с напряженностью

Ео ~ 3 х 104 В/см. Моделирование действия полей с еще большими значениями не представляет принципиальных трудностей и требует только уверенности в адекватности используемой модели реальной физической системе.

Переход же к более низким значениям напряженности электрического поля ведет к усложнению задачи достаточно точного воспроизведения /^ поскольку область её локализации быстро уменьшается. И даже в присутствии поля в условиях использования арифметики двойной точности область доступных для воспроизведения значений функции распределения оказывается ограниченной. Для оценки нижней границы доступных для моделирования полей был выполнен ряд численных экспериментов в диапазоне значений 3 х 104 < Ео < 3 В/см. Результаты для /(р1,р2,^ = 0) при Е0 = 3 х 102 В/см и Е0 = 3 В/см показаны на рисунках 5 и 6.

1(РьР2

Рисунок 5. Функция распределения f (р1,р2,1 = 0) для момента времени с максимальным значением напряжённости внешнего поля 3 х 102 В/см. Остальные параметры аналогичны рисунку 3

Использовались точно такие же сетки, как и при получении результатов, представленных на рисунке 3.

¡(РрР^

Рисунок 6. Функция распределения f (р1,р2,Ь = 0) для момента времени с максимальным значением напряжённости внешнего поля 3 В/см. Остальные параметры аналогичны рисункам 3 и 5

Уже при Е0 = 3 х 102 В/см и в присутствии внешнего поля в значительной части области —п < р1 < п, —п < р2 < п значения функции распределения оказываются меньше порога 1.0 х 10-14. Но доступная для воспроизведения часть значений функции распределения ведёт себя гладко, достигая максимальных значений ~ 1.25 х 10-7 в симметрично расположенных точках р1 ± 0.025, р ± 0.025. Это дает основания полагать, что определяемый некоторым пороговым значением функции распределения контур будет охватывать всю область локализации /(Р1,Р2, ^ ^ то). При выборе в качестве порогового значения /(р1,р2, 0) = 1.0 х 10-7 локализуется область —0.05 < р1 < 0.05, —0.05 < р2 < 0.05 за пределами которой значений /(р1 ,р2, 0) > 1.0 х 10-7 не существует. Такая оценка позволяет уменьшить площадь исследуемой области импульсного пространства с использованием более плотных сеток примерно в 4000 раз. Конкретный выбор используемого порогового значения, его соотнесение с границей доступных для вычисления значений /(р1,р2,^ ^ то), зависимость этого порогового

значения от характера зависимости электрического поля от времени требуют дальнейших исследований.

На нижней границе Е0 = 3 В/см рассматривавшегося диапазона область доступных для вычисления значений очень мала, а максимальные значения функции распределения составляют всего /(р1,р2, 0) = 1.255 х 10-11, которые достигаются в тех же точках Р1 ± 0.025, р2 ± 0.025 (в использовавшейся сетке нет точек, расположенных ближе к началу координат). Тем не менее, этого достаточно для локализации функции распределения и точного её воспроизведения.

Представленные результаты показывают, что во всем рассматривавшемся диапазоне 3 < Ео < 3 х 104 В/см обеспечивается возможность использования результатов вычисления функции распределения в присутствии внешнего поля для определения области локализации её максимальных значений.

3. Воспроизведение динамики наблюдаемых характеристик

Процедура моделирования позволяет получит исчерпывающую информацию как о /(рьр>2,^ ^ то), так и о значениях функции распределения в любой промежуточный момент времени /(р1,р2,^) начиная с момента включения внешнего поля. В общем случае это достаточно объёмные трехмерные массивы, данные из которых необходимо затем анализировать. При этом содержащаяся в /(р1,р2,^) для каждого конкретного момента времени информация избыточна, так как в полном объеме не верифицируема экспериментально. Для оценки наблюдаемых параметров моделируемых процессов достаточно интегральных характеристик этой функции вида:

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

(2)

(3)

зависит от времени. Интегрирование выполняется по области определения f (pi,p2), а фактически в пределах области, определяющей основной вклад в наблюдаемые интегральные характеристики.

Если мы точно знаем набор интегральных характеристик, поведение которых нас интересует, разделение процедуры вычисления и сохранения 3D массива значений f (pi,p2,t) и процедуры интегрирования оказывается не целесообразным. В версию моделирующей системы, реализованную с использованием MPI, был введен вариант работы, обеспечивающий вычисление (и сохранение) для последовательности временных точек только заданного набора интегральных значений. Для этого используется процедура редукции, реализуемой средствами MPI. Это обеспечивает очень существенное сокращение требуемых ресурсов на хранение и последующую обработку данных, поскольку вместо 3D массива мы имеем дело только с несколькими 1D массивами задаваемой длины.

В качестве примера на рисунке 7 приведены результаты вычисления поверхностной плотности носителей зарядов (2) и поверхностной плотности тока (3) в условиях действия импульса внешнего поля с максимальным значением напряженности 3 х 102 В/см. Остальные параметры те-же, что использовались ранее. Эти массивы содержат по 361 значению. Каждое из значений соответствует определенному моменту времени с шагом примерно 5 фемтосекунд. На рисунках представлено поведение наблюдаемых параметров с начала появления внешнего поля (когда носители и ток отсутствуют) и до их выхода на фиксированные значения после выключения внешнего поля.

300 225

150 _ g

75 > 0 BÏ -75 -150

300 225 150 75 0

-75 -150

t [s]

Рисунок 7. Поверхностная плотность носителей зарядов п(Ь) (2) и поверхностная плотность тока j(t) (3) в условиях действия импульса внешнего поля с максимальным значением напряженности 3 х 102 В/см

2x10

1.5x10

1x10

5x10

0

-5x10

-1x10

SIMULATE THE BEHAVIOR OF GRAPHENE IN EXTERNAL ELECTRIC FIELDS

Текст научной работы на тему «МОДЕЛИРОВАНИЕ ПОВЕДЕНИЯ ГРАФЕНА ВО ВНЕШНИХ ЭЛЕКТРИЧЕСКИХ ПОЛЯХ»