CC BY
4
.RfcUjN ПРОГРАММНЫЕ СИСТЕМЫ: ТЕОРИЯ И ПРИЛОЖЕНИЯ ISSN 2079-3316
научная статья математическое моделирование
УДК 519.688:519.876.5
10.25209/2079-3316-2023-14-2-27-47
Моделирование поведения двухуровневой квантовой системы с использованием масштабируемых регулярных сеток
Анатолий Дмитриевич Панферов1", Наталия Владимировна Поснова2, Анастасия Алексеевна Ульянова3
Саратовский государственный университет им. Н.Г. Чернышевского, Саратов, Россия
Аннотация. Представлены результаты программной реализации покрытия двумерной области допустимых состояний моделируемой двухуровневой квантовой системы регулярными сетками различного масштаба. Метод был предложен и реализован применительно к задаче воспроизведения поведения графена во внешних электрических полях с использованинем квантового кинетического уравнения. Численное решение кинетического уравнения позволяет точно и полно воспроизводить поведение модели, однако проблемой является минимизация необходимых вычислительных ресурсов при обеспечении требуемого качества вычисления наблюдаемых параметров. Сложность задачи прямо пропорциональна количеству узлов сетки в двумерном пространстве состояний, в каждом из которых для воспроизведения эволюции заселенности верхнего энергетического уровня решается система обыкновенных дифференциальных уравнений. Исследование особенностей поведения модели позволило разработать и реализовать экономичную процедуру построения рассчетной сетки, обеспечивающую локализацию области формирования возбуждённых состояний и точное воспроизведение наблюдаемых параметров процесса.
Ключевые слова и фразы: численное моделирование, графен, функция распределения носителей заряда, оптимальный выбор расчетной сетки
Благодарности: Исследование выполнено за, счет гранта. Российского научного фонда, № 23-21-00047mL
Для цитирования: Панферов А. Д., Поснова Н.В., Ульянова А. А. Моделирование поведения двухуровневой квантовой системы с использованием масштабируемых регулярных сеток // Программные системы: теория и приложения. 2023. Т. 14. № 2(57). С. 27-47. https://psta.psiras.ru/read/psta2023_2_27-47.pdf
© Панферов А. Д., Поснова Н.В., Ульянова А. А.
2023
Введение
Активный поиск новых материалов для современной электроники повлек за собой развитие исследований структур с пониженной размерностью. Технологически удобными объектами такого типа являются двумерные пленки, способные устойчиво существовать в свободном пространстве или на подложках. Свойства электронов в них и особенности процессов их взаимодействия с внешними электрическими полями могут быть описаны только с привлечением принципов и подходов квантовой механики. Наиболее известным примером материалов этого типа является графен, доступ к образцам которого позволил разработать и экспериментально верифицировать уже ставшую классической модель безмассовых фермионов [1,2]. В настоящее время интерес к этому материалу обусловлен особенностями его отклика на действие внешних электрических полей [3-8].
Для описания процессов образования и эволюции носителей заряда в двумерном монослое графена в условиях действия внешнего классического электрического поля в работах [9-11] был предложен и развит подход, основанный на использовании квантового кинетического уравнения. Метод квантового кинетического уравнения достаточно универсален и может использоваться при рассмотрении других аналогичных материалов. В данном подходе подсистема слабосвязанных электронов материала рассматривается как квантовая система с двумя энергетическими уровнями. Конкретный вид закона дисперсии (зависимости энергии состояний от импульса) может быть различным. В простейшем случае он соответствует приближению безмассовых фермионов, но может использоваться и строгая модель взаимодействия ближайших соседей. Кинетическое уравнение определяет эволюцию функции распределения f (р,£) носителей по возбуждённым состояниям в условиях внешнего воздействия и может быть представлено в виде системы обыкновенных дифференциальных уравнений, определяемой через гамильтониан модели.
Трудности использования этого подхода обусловлены тем обстоятельством, что лежащая в его основе система уравнений может решаться только численными методами и для каждой рассматриваемой точки импульсного пространства решение необходимо искать независимо. Это приводит к высокой вычислительной сложности процедуры получения наблюдаемых характеристик. В работе [12] была продемонстрирована возможность использования адаптивной сетки на основе квадродерева для вычисления конечного состояния /(р, £ ^ то), формирующегося
после завершения действия внешнего электрического поля с импульсным характером зависимости от времени.
Использование сетки переменной плотности в области локализации максимальных значений функции распределения обеспечивает существенную экономию ресурсов по сравнению с использованием регулярных сеток. Однако проблема предварительного определения местоположения области локализации максимальных значений не имела эффективного решения. На основе практики моделирования процессов в рассматриваемых системах было обнаружено, что локализацию областей концентрации возбуждённых состояний (где требуется использовать плотные сетки) удобнее выполнять по данным о / (р, £) не после завершения внешнего воздействия, а во время его. Это было продемонстрировано в работе [13], в которой сопоставлялись особенности поведения в импульсном пространстве /(р, £ ^ то) и / (р, £Етах), где tEmax определялось как момент времени с максимальным значением внешнего воздействия (напряжённости внешнего электрического поля).
В развитии этой идеи для моделей рассматриваемого типа в представляемой работе выполнены исследования и сопоставление поведения в импульсном пространстве функции распределения в конечном состоянии /(р, £ ^ то), в момент максимума внешнего воздействия /(р,^Етах), а также максимальных значений функции распределения в выбранной точке импульсного пространства за все время действия внешнего поля /тах(р) и средних значений функции распределения за время действия внешнего поля /ауд (р). Это позволило разработать итерационный метод построения разномасштабных регулярных сеток, обеспечивающий в рамках унифицированной процедуры локализацию области генерации возбуждённых состояний с охватом всех рассматриваемых значений (р) (полного образа первой зоны Бриллюэна) и формирование в ней плотных сеток разного масштаба.
1. Модель и её параметры
При исследовании особенностей реакции перспективных двумерных материалов на действие внешних электрических полей эффективно одноэлектронное приближение. В этом приближении квантовая эволюция электронов, не связанных жестко с конкретными атомами и имеющими возможность реагировать на воздействие внешнего поля, определяется уравнением Шрёдингера через одночастичный оператор Гамильтона, отражающий специфику внутрнних полей кристаллической решетки.
В отсутствии явной зависимости гамильтониана от времени в импульсном представлении уравнение Шрёдингера принимает стационарную форму
(1) нрФ(р,г) = е(р)Ф(р,г),
определяющую спектр состояний (зонную структуру) е(р) и соответствующие этим состояниям собственные функции Ф(р,Ь). Эволюция стационарных состояний во времени определяется в явной форме фазовым множителем ехр(-ге(р)Ь/К).
Появление классического внешнего электрического поля выводит систему из стационарного состояния и меняет характер её эволюции. Действие классического поля на квантовую систему описывается через модификацию гамильтониана путем замены в нем оператора импульса на квазиимпульс по правилу
(2) Рк ^ Рк = Рк - еАк(Ь), (к = 1,22),
где Ак (Ь) — компоненты векторного потенциала внешнего поля и мы принимаем, что заряд электронов отрицателен, а его абсолютное значение совпадает с величиной элементарного заряда е. Наблюдаемая напряжённость электрического поля Ек (Ь) определяется выражением
Ек (Ь) = -А к (Ь).
Воспроизведение эволюции в условиях явной зависимости гамильтониана от времени для реальных систем является вычислительно сложной задачей, решаемой только численными методами. При этом удобнее работать не с самой одночастичной волновой функцией, а с определяемой через неё функцией распределения ](р, Ь) электронов по одночастичным энергетическим состояниям. Её эволюция может быть определена через кинетическое уравнение вида:
0(Р,ь) = 2^(Р,ь)п(р,г),
(3) I й(р,ь) = \(р,г)(1 - 2/(р,ь)) - у(р,г), [цр,г) = ЩА й(р,ь),
где й(р,Ь) и у(р,Ь) вспомогательные функции (детали можно посмотреть в работе [11]). Поведение моделируемой системы определяется переменными коэффициентами е(р,Ь) и Х(р,Ь). Первый из них выражается через стационарный спектр е(р) из уравнения (1) путем замены (2). Второй также определяется через исходный гамильтониан системы с использованием (2).
Так, при рассмотрении двумерной гексагональной решетки графена, представляющей суперпозицию двух одинаковых подрешеток, общий вид
гамильтониана имеет вид
(4) Н(„)= ( ^ )
с двумя собственными значениями л/БВ и -л/ВВ*. Соответственно, в этом случае е(р, £) = \/В(р,1)В*(р,1), а выражение для определения второго коэффициента из системы уравнений (3) имеет вид
(5) = В'(р,г)В"(р,г) - Б'(р,г)В"(р,г)
(5) ^г) = ЩР^Т) ,
где В'(р,£) и В''(р,£) — вещественная и мнимая части матричного элемента гамильтониана (4):
(6) В(р,г) = В'(р,г) + 1Вп(р,г).
В реальных случаях, например, при использовании модели взаимодействия ближайших соседей, явный вид этих выражений может быть достаточно сложным.
Периодической пространственной структуре моделируемого материала соответствует периодическая обратная решетка. Для случая графена её вид представлен на рисунке 1:
Рисунок 1. Структура обратной решетки графена. Единица измерения [р]=—^ где ^ = 0.142 нм—длина ребра шестигранной решетки в прямом пространстве.
При моделировании необходимо воспроизвести эволюцию функции распределения (решить систему уравнений (3) при всех рассматриваемых значениях p) в первой зоне Бриллюэны (желтый шестиугольник) или (в силу периодичности) эквиваленте этой зоны. В качестве эквивалента была выбрана элементарная ячейка обратной решетки, построенная на векторах трансляций Ъ\ и 62 с центром в точке M, которая будет использоваться в качестве начала системы координат. При таком выборе области определения для функции распределения в ней оказываются две неэквивалентные точки Дирака K и K', расположенные симметрично относительно начала отсчета системы координат на оси ординат.
Для реализации возможностей модели разработан и развивается комплекс программ. Отработка его алгоритмов и функционала осуществляется поэтапно. Вначале решения тестируются на макете системы, реализованном средствами Wolfram Mathematica. Затем оформляются в виде модулей на C/C++.
Наиболее ресурсоемким является процесс решения кинетического уравнения (3), представляющего из себя систему обыкновенных дифференциальных уравнений. На этапе работы с макетом для этого используются хорошо зарекомендовавшие себя решатели Mathematica. Модули на C реализованы с использованием библиотеки GSL.
Разработка вычислительного модуля, ориентированного на использование ресурсов графических сопроцессоров NVIDIA, потребовала перехода на решатели Odeint из библиотеки Boost C+—Ь, хорошо интегрируемые в программы на CUDA. Все эти решатели могут использовать различные алгоритмы. Проведённое тестирование показало, что в случае рассматриваемого уравнения преимущество имеют явные методы высокого порядка, основанные на классе методов Рунге-Кутты: runge_kutta_cash_karp54, runge_kutta_dopri5, runge_kutta_fehlberg78.
Важным с точки зрения организации работы программы является возможность решать системы уравнений (3) в разных точках p независимо (параллельно). Такой параллелизм по данным без обмена промежуточными результатами позволяет эффективно использовать ресурсы современных многоядерных процессоров, графических ускорителей и мультикомпьюте-ров. На уровне макета используются возможности распараллеливания самой Mathematica. Решения на C/C++ используют MPI, что позволяет запускать их на мультикомпьютерах.
Наблюдаемые физические параметры воспроизведение которых является целью моделирования, определяются через усреднение соответствующих величин вида
по вычисленной функции распределения, где интегрирование необходимо выполнить в пределах элементарной ячейки обратной решетки. В отсутствии аналитического представления / (р,1) это интегррование может быть выполнено только численно. Но, если зависимость /(р,£) от времени при заданном значении р определена системой уравнений (3), информации о поведении функции распределения в двумерном импульсном пространстве при фиксированном значении £ модель в явной форме не содержит. По этой причине возникает задача о покрытии области определения /(р,Ь) в импульсном пространстве некоторой сеткой, определив на которой значения функции распределения путем решения системы уравнений (3) в каждом узле, можно получить достаточно точные оценки для наблюдаемых (7). Предложенный вариант решения этой задачи излагается далее.
2. Локализация области генерации возбуждённых состояний по характеристикам функции распределения носителей
Максимальная энергия возбуждённых состояний достигает в графене примерно 9 эВ. При этом в окрестностях точек Дирака она близка к нулю. Это приводит к крайне неоднородному распределению возбуждённых состояний. Подробно эта проблема обсуждалась в работе [13]. Там же была продемонстрирована принципиальная разница в поведении / (р) в конечном состоянии после выключения поля и в его присутствии.
На основании отмеченной особенности было принято решение исследовать характер поведения в импульсном пространстве таких характеристик функции распределения, как среднее значение за время возмущающего воздействия и максимальное её значение за время действия возмущающего воздействия. Обозначим их /ащд(р) и /тах(р). На рисунке 2 показано в высоком разрешении поведение функции распределения в окрестностях одной из точек Дирака для короткого импульса электромагнитного поля с круговой поляризацией. Области, в которых присутствуют достаточно большие значения из интервала [0.1 ^ /(р) ^ 1.0], совпадают.
(7)
Рисунок 2. Функция распределения после завершения действия возмущения f (р, £ ^ то), её усреднённые за время действия возмущения значения favg (р) и максимальные зарегистрированные значения /тах(р) в окрестностях нижней точки Дирака К (см. рисунок 1).
Сравним поведение представленных функций в глобальном масштабе. На рисунке 3 показаны /ауд (р) и ¡.тах (р) для отрицательных значений р1. В силу очень большого диапазона значений выполнено приведение к логарифмическому масштабу преобразованием вида / ^ 10 /.
Ожидаемо /ауд(р) < /тах(р), но при этом их поведение во всей рассматриваемой области гладкое и очень похожее /ауд(р) ~ /тах(р). Исходя из этого подобия далее будем использовать ¡тах(р) и сравним её
Рисунок 3. Глобальное поведение значений /аьд(р) и /тах(р).
с функцией распределения для конечного состояния модели /(р,Ь ^ то). Такое сравнение приведено на рисунке 4 с аналогичным приведением к логарифмическому масштабу.
Рисунок 4. Глобальное поведение значений /тах(р) и функции распределения для конечного состояния модели /(р^ <).
Представленные результаты полностью согласуются с выводами,
сделанными ранее в работе [13]. Корректное воспроизведение f (p,t ^ ж) возможно только в ограниченных областях импульсного пространства. Вне этих областей значения функции недоступны для вычисления при использовании 64-битной арифметики. При этом имеет место строгое совпадение мест локализации максимумов f (p, t ^ ж) и fmax(p), что позволяет использовать последнюю в качестве удобного инструмента для определения положения максимумов.
При ограничении условием fmax(p) < 0.1 эта функция ведет себя гладко во всей области определения. Для локализации по её значениям областей с максимумами будем использовать регулярную сетку N х N, определяющую в импульсном пространстве дискретный набор точек [pi,pj}. На сетке будут отбираться узлы с fmax(Pi,Pj) > fiim, где fHm — нижняя граница интересующего нас диапазона значений. Уже на регулярной сетке с N = 32 (Api = Лр2 = 4п/32) для использовавшихся выше параметров импульса локализация уверенно выполняется при f¡im = 1.0 х 10-9. На рисунке 5 представлено распределение значений на такой сетке для отрицательного полупространства pi ^ 0.
Рисунок 5. Распределение значений 10log(/max(pi,pj)) при N = 32.
На рисунке 6 выделены только значения, удовлетворяющие условию /тах(рг,рз) > 1.0 х 10-9 для всей области допустимых значений (р1,р2).
Рисунок 6. Выборка значений 10\og(fmax(pi,р^)) для /тах (Р1,РЗ ) > 1.0 X 10-9 при N = 32.
Необходимо отметить, что выполненная процедура не предоставляет точных координат области, показанной на рисунке 2. Пусть критерий попадания в эту область 0.01 < /тах(рг,рз) < 1.0. Максимальное же значение /тах(рг,рз) из полученных на сетке 32 х 32 и показанных на рисунках 5 и 6 всего 3.99 х 10-7. Возможное использование на первом этапе сеток с меньшим шагом принципиально не решает эту проблему. Представленные результаты получены только для одного из многих возможных сочетаний параметров моделируемого процесса, но они вполне показательны.
3. Выборочное масштабирование сетки
Продемонстрированная выше процедура позволяет уменьшить область поиска до нескольких ячеек исходной сетки. Конкретное их число будет
определяться выбранным уровнем отсечения /цт. При /цт = 1.0 х 10-9 значение fmax(pi,pj) превышает этот уровень в 34 узлах сетки. При /ит = 1.0 х 10-7 пороговое значение будет превышено только в 4 узлах. Эти узлы расположены попарно в двух разнесенных областях.
Покажем, что этого достаточно для локализации области, где /тах(р) имеет близкие к единице значения. Дальнейшее исследование поведения /тах(р) в окрестностях выбранных точек выполним по аналогии с первым шагом, разбивая каждую ячейку с соответствующей узловой точкой в середине регулярной сеткой М х М. Сохраним значение М = 32. Результаты показаны на рисунке 7. Четыре ячейки исходной сетки (далее будем говорить о ней как о сетке 0-го уровня) объединены попарно.
Рисунок 7. Распределение значений 10log(/max(p;,Pj)) для двух пар ячеек нулевого уровня.
На рисунке 8 представлена общая картина, получающаяся при объединении новых результатов с результатами для узлов нулевого уровня:
На показанных на рисунке 7 сетках первого уровня присутствует 32 узла со значениями fmax(pi,pj) ^ 0.01. Продолжим процедуру и выполним на них следующий шаг дискретизации с использованием
Рисунок 8. Совмещение результатов, полученных для /тах (Рг,Рз) на сетках нулевого и первого уровней.
масштабированной регулярной сетки 32 х 32. На рисунке 9 показаны результаты для ячейки первого уровня, соответствующей максимуму на правой части Рисунка 7.
Рисунок 9. Распределение /тах(рг,р]) в ячейке сетки первого уровня.
На этом этапе становится доступна мелкомасштабная структура
функции распределения, причем все представленные значения лежат в диапазоне 0.1 < fmax(Pi,Pj) < 1.0.
Полная картина для одного из максимумов, показанных на рисунках 7 и 8, сгруппирована из данных для 16 ячеек сетки первого уровня на рисунке 10. На нем представлены значения 0.01 < fmax(pi,pj) < 1.0,
N = 32.
приведённые к логарифмическому масштабу.
В результате двух итераций получена гибридная сетка с тремя масштабами покрытия исследуемой области с отношением размеров ячеек 1 : 1/32 : 1/32/32. Общее количество узлов гибридной сетки ~ 3.7 х 104. Процедура обеспечивает как локализацию областей с максимальными значениями функции распределения, определяющими наблюдаемые параметры моделируемого процесса, так и точное воспроизведение мелкомасштабной структуры функции распределения в двух неэквивалентных областях.
Располагая такой сеткой, можно экономично и точно исследовать поведение модели на всем протяжении процесса её взаимодействия с импульсом внешнего поля. На рисунке 11 показаны функции распределения для ряда последовательных моментов времени {¿о,Ь\,Ь2, ¿з}. Первый из них соответствует максимуму напряженности внешнего электрического
■
I
Рисунок 10. Распределение значений 10log(/max(p;,pj)) при
Рисунок 11. Воспроизведение поведения /(р ,р2,£) на примере четырех последовательных моментов времени.
поля, последний — полному выключению внешнего воздействия, {¿1,^2 } — промежуточные значения. Приведены результаты для области покрытия сеткой второго уровня, показанной на рисунке 10. Во всех случаях обеспечивается воспроизведение тонкой структуры распределения возбуждённых состояний и полный охват области их локализации.
Оценку эффективности предложенной процедуры можно сделать сравнивая построенную сетку переменного масштаба с эквивалентной ей по обеспечиваемой детализации равномерной сеткой. В области воспроизведения тонкой структуры функции распределения в окрестностях её максимумов (например, рисунок 9) расстояние между узлами построенной сетки составляет 4п/32/32/32. Если воспользоваться равномерной сеткой с таким шагом, для покрытия всей области определения функции распределения в импульсном пространстве потребуется ~ 3.1 х 108 узлов, что почти на четыре порядка больше числа узлов в построенной
сетке переменного масштаба. Обеспечиваемая экономия вычислительных ресурсов делает возможным моделирование процессов рассматриваемого типа на современных высокопроизводительных рабочих станциях.
Фактические затраты времени на каждый узел сетки зависят от параметров моделируемого процесса, программной реализации, используемой аппаратной платформы и ряда других факторов. Демонстрация процедуры построения сетки выполнена при фиксированных параметрах моделируемого процесса на системе с CPU Intel(R) Core(TM) Í5-11400 с б физическими ядрами и базовой тактовой частотой 2.60 GHz, Кэш L3 12.0 MB, объём памяти З2.0 GB. Все шаги процедуры, включая решение системы уравнений (3), были реализованы с использованием возможностей Mathematica. Характерные затраты времени на воспроизведение поведения модели в каждом узле строящейся сетки оцениваются в O.l239 c при вычислениях в последовательном режиме (используется одно ядро) или O.O2O65 c при параллельной работе на б физических ядрах. Соответственно, оценка времени на воспроизведение поведения модели с использованием всей построенной сетки составляет 764. O c при полном использовании ресурсов системы.
Заключение
Численное моделирование многочастичных систем, поведение которых в определяющей мере подчиняется законам квантовой механики, ресурсоемко и требует оптимизации всех элементов процедуры. В представленной работе продемонстрирована реализация метода построения расчетной сетки для покрытия области определения функции распределения носителей заряда в двумерном импульсном пространстве модельной системы с двумя энергетическими уровнями. Такая модель позволяет описывать поведение перспективных двумерных материалов в условиях действия на них внешних электрических полей. Например, графена, причем как в наиболее распространенном упрощении приближения безмассовых фермионов, так и в более строгой модели взаимодействующих ближайших соседей. Предложенный метод позволяет в рамках одной процедуры совместить поиск области локализации возбуждений в импульсном пространстве модели и воспроизведение детальной картины их распределения в импульсном пространстве и во времени.
Список литературы
[1] Novoselov K. S., Geim A.K., Morozov S.V., Jiang D., Katsnelson M. I.,
Grigorieva, I. V., Dubonos S. V., Firsov A. A. Two-dimensional gas of massless dirac fermions in graphene // Nature - 2005 - Vol. 438.- Pp. 197-200. 28
[2] Castro Neto A. H., Guinea F., Peres N. M. R., Novoselov K. S„ Geim A. K. The eletronie properties of graphene // Rev. Mod. Phys.- 2009.- Vol. 81.- No. 1.- 109.
d 28
[3] Glazov M. M., Ganichev S. D. High frequency electric field induced nonlinear effects in graphene jj Physics Reports.- 2014,- Vol. 535,- No. 3,- Pp. 101-138. 28
[4] Bowlan P., Martinez-Moreno E., Reimann K., Elsaesser T., Woerner M. Ultrafast terahertz response of multilayer graphene in the nonperturbative regime // Phys. Rev. В.- 2014,- Vol. 89,- 041408. 28
[5] Baudisch M., Marini A., Cox J. D., Zhu Т., Silva F., Teichmann S., Massicotte M., Koppens F., Levitov L. S., Garcia de Abajo F. J., Biegert J. Ultrafast nonlinear optical response of Dirac fermions in graphene // Nature Communications.- 2018.-Vol. 9,- 1018,- 6 Pp. 28
[6] Zi-Yu Chen, Rui Gin Circularly polarized extreme ultraviolet high harmonic generation in graphene // Optics Express.- 2019.- Vol. 27.- No. 3.- Pp. 3761-3770. d 28
[7] Wenwen Mao, Angel Rubio, Shunsuke A. Sato Terahertz-induced high-order harmonic generation and nonlinear charge transport in graphene // Phys. Rev. B.-2022,- Vol. 106,- 024313. 28
[8] Soonyoung Cha, Minjeong Kim, Youngjae Kim, Shinyoung Choi, Sejong Kang, Hoon Kim, Sangho Yoon, Gunho Moon, Taeho Kim, Ye Won Lee, Gil Young Cho, Moon Jeong Park, Cheol-Joo Kim, Kim B. J., JaeDong Lee, Moon-Ho Jo, Jonghwan Kim Gate-tunable quantum pathways of high harmonic generation in graphene // Nature Communications.- 2022,- Vol. 13,- 6630,- 10 Pp. 28
[9] Smolyansky S. A., Churochkin D. V., Dmitriev V. V., Panferov A. D., Kampfer B. Residual currents generated from vacuum by an electric field pulse in 2+1 dimensional QED models // EPJ Web of Conferences.- 2017,- Vol. 138,- 06004,5 Pp. 28
[10] Panferov A. D., Smolyansky S. A., Blaschke D. B., N. Gevorgyan T. Comparing two different descriptions of the I-V characteristic of graphene: theory and experiment // EPJ Web of Conferences.- 2019,- Vol. 204,- 06008,- 6 Pp. I
[11] Smolyansky S. A., Panferov A. D., Blaschke D. B., Gevorgyan N. T. Nonperturbative kinetic: description of electron-hole excitations in graphene in a time dependent electric field of arbitrary polarization // Particles.- 2019.- Vol. 2.- No. 2.-Pp. 208-230. 28 so
[12] Панферов А. Д., Маханьков А. В., Трунов А. А. Использование адаптивной сетки на основе квадродерева для моделирования конечного состояния квантово-полевой системы при импульсном внешнем воздействии // Программные системы: теория и приложения.- 2020.- Т. 11.- № 1.- С. 79-92 .
, U R L % 28
[13] Панферов А. Д., Новиков Н. А., Трунов А. А. Моделирование поведения графена во внешних электрических полях // Программные системы: теория и приложения.- 2021.- Т. 12.- С. 3. url 2в зз 36
Поступила в редакцию одобрена после рецензирования принята к публикации опубликована онлайн
20.02.2023; 12.05.2023; 04.06.2023; 07.07.2023.
д.т.н. А. М. Цирлин
Рекомендовал к публикации
Информация об авторах:
Анатолий Дмитриевич Панферов К.ф.-м.н., зам. начальника УЦИТ Саратовского государственного университета им. Н. Г. Чернышевского по научно-производственной деятельности. Научные интересы: высокопроизводительные вычисления, параллельное программирование, численное решение квантовых кинетических уравнений, моделирование процессов вакуумного рождение частиц в КЭД, генерации носителей в полупроводниках в том числе бесщелевых, процессов на ранних стадиях столкновения релятивистских ядер
e-mail:
0000-0003-2332-0982 [email protected]
Наталия Владимировна Поснова
Магистрант Саратовского государственного университета им. Н. Г. Чернышевского. Научные интересы: моделирование физических процессов на высокопроизводительных вычислительных системах, параллельное программирование.
¡ГО 0000-0003-1259-1867 e-mail: [email protected]
Анастасия Алексеевна Ульянова
Магистрант Саратовского государственного университета им. Н. Г. Чернышевского. Научные интересы: моделирование физических процессов на высокопроизводительных вычислительных системах, параллельное программирование.
e-mail:
0000-0001-9519-9822 [email protected]
Вклад авторов: А.Д. Панферов — 99.95% (идея, методология, создание стилевого файла, сбор материала, написание черновой версии, доработка и редактирование, визуализация, администрирование); Н.В. Поснова — 0.05% (доработка и редактирование); А. А. Ульянова — 0.05% (доработка и редактирование).
Авторы заявляют об отсутствии конфликта интересов.
BUS PROGRAM SYSTEMS: THEORY AND APPLICATIONS ISSN 2079-3316
Research Article mathematical modeling
UDC 519.688:519.876.5
10.25209/2079-3316-2023-14-2-27-47
Simulation the behavior of a two-level quantum system using scalable regular grids
Anatolii Dmitrievich Panferov1, Natalia Vladimirovna Posnova2, Anastasia Alekseevna Ulyanova3
Saratov State University, Saratov, Russia
Abstract. The paper demonstrates a method for constructing a grid with a changing scale for modeling the behavior of a two-level quantum system. The method was proposed and implemented in relation to the problem of reproducing the behavior of graphene in external electric fields using the quantum kinetic equation. A direct numerical solution of the kinetic equation makes it possible to accurately and completely reproduce the behavior of models of this type. But a large amount of computation is required for a good quality reproduction of the observed parameters. The complexity of the problem is directly proportional to the number of grid nodes in the two-dimensional space of states, in each of which a system of ordinary differential equations is solved to reproduce the evolution of the population of the upper energy level. The study of the behavior of the model made it possible to develop and implement an economical procedure for constructing a computational grid, which ensures the localization of the region of formation of excited states and exact reproduction of the observed process parameters. (In Russian).
Key words and phrases: numerical simulation, graphene, charge carrier distribution function, optimal choice of computational grid
2020 Mathematics Subject Classification: 65M50; 81T10, 81T40
Acknowledgments: The study was supported by a grant from the Russian Science Foundation No. 23-21-00047"
For citation: Anatolii D. Panferov, Natalia V. Posnova, Anastasia A. Ulyanova. Simulation the behavior of a two-level quantum system using scalable regular-grids . Program Systems: Theory and Applications, 2023, 14:2(57), pp. 27-47. (In Russ.). https://psta.psiras.ru/read/psta2023_2_27-47.pdf
© Panferov A. D., Posnova N. V., Ulyanova A. A.
2023
References
[1] K. S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, M.I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, A. A. Firsov. "Two-dimensional gas of massless dirac fermions in graphene", Nature, 438 (2005), pp. 197-200.
[2] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, A. K. Geim. "The eletronic properties of graphene", Rev. Mod. Phys., 81:1 (2009), id. 109.
[3] M. M. Glazov, S. D. Ganichev. "High frequency electric field induced nonlinear effects in graphene", Physics Reports, 535:3 (2014), pp. 101-138.
[4] P. Bowlan, E. Martinez-Moreno, K. Reimann, T. Elsaesser, M. Woerner. "Ultrafast terahertz response of multilayer graphene in the nonperturbative regime", Phys. Rev. B, 89 (2014), id. 041408.
[5] M. Baudisch, A. Marini, J.D. Cox, T. Zhu, F. Silva, S. Teichmann, M. Massicotte, F. Koppens, L. S. Levitov, de Abajo F. J. Garcia, J. Biegert. "Ultrafast nonlinear optical response of Dirac fermions in graphene", Nature Communications, 9 (2018), id. 1018, 6 pp.
[6] Chen Zi-Yu, Gin Rui. "Circularly polarized extreme ultraviolet high harmonic generation in graphene", Optics Express, 27:3 (2019), pp. 3761-3770.
[7] Mao Wenwen, Rubio Angel, A. Sato Shunsuke. "Terahertz-induced high-order harmonic generation and nonlinear charge transport in graphene", Phys. Rev. B, 106 (2022), id. 024313.
[8] Cha Soonyoung, Kim Minjeong, Kim Youngjae, Choi Shinyoung, Kang Sejong, Kim Hoon, Yoon Sangho, Moon Gunho, Kim Taeho, Won Lee Ye, Young Cho Gil, Jeong Park Moon, Kim Cheol-Joo, B.J. Kim, Lee JaeDong, Jo Moon-Ho, Kim Jonghwan. "Gate-tunable quantum pathways of high harmonic generation in graphene", Nature Communications, 13 (2022), id. 6630, 10 pp.
[9] S.A. Smolyansky, D.V. Churochkin, V.V. Dmitriev, A.D. Panferov, B. Kampfer. "Residual currents generated from vacuum by an electric field pulse in 2+1 dimensional QED models", EPJ Web of Conferences, 138 (2017), id. 06004, 5 pp.
[10] A.D. Panferov, S.A. Smolyansky, D.B. Blaschke, N. T. Gevorgyan. "Comparing two different descriptions of the I-V characteristic of graphene: theory and experiment", EPJ Web of Conferences, 204 (2019), id. 06008, 6 pp. d
[11] S.A. Smolyansky, A.D. Panferov, D.B. Blaschke, N.T. Gevorgyan. "Nonperturbative kinetic description of electron-hole excitations in graphene in a time dependent electric field of arbitrary polarization", Particles, 2:2 (2019), pp. 208-230.
[12] A. D. Panferov, A. V. Makhankov, A. A. Trunov. "The use of an adaptive mesh based on a quadtree for modeling the final state of a quantum field system under pulsed external action", Program Systems: Theory and Applications, 11:1 (2020), pp. 93-105. URL
[13] A. D. Panferov, N. A. Novikov, A. A. Trunov. "Simulate the behavior of graphene in external electric fields", Program Systems: Theory and Applications, 12 (2021), pp. 3 (in Russian), .url
