CC BY
12
ББК 32.973
УДК 519.688, 004.942
А. Д. Панферов, А. В. Маханьков
Моделирование действия коротких оптических импульсов на графен
Аннотлция. Особенности взаимодействия высокочастотных импульсных электрических полей с графеном в настоящее время являются предметом интенсивных исследований. В работе представлены результаты апробации программной системы моделирования таких процессов на примере ультракоротких лазерных импульсов оптического диапазона с различной поляризацией. Разрабатываемая авторами программная система построена на базе нового теоретического подхода с использованием квантового кинетического уравнения. Ключевым элементом модели является система обыкновенных дифференциальных уравнений с нелинейно зависящими от времени и параметров задачи коэффициентами.
Необходимость анализировать поведение решений этой системы уравнений в области изменения нескольких параметров ведёт к полиномиальной сложности вычислений. Неизвестность характера параметрической зависимости решений требует нескольких итераций выбора покрывающих сеток. Система моделирования адаптирована для использования на массово параллельных вычислительных системах.
Ключевые слова и фразы: численное моделирование, высокопроизводительные вычисления,
графен, кинетическое уравнение.
Введение
В настоящее время интенсивно исследуются особенности взаимодействия оптического излучения с новыми двумерными или псевдо-двумерными материалами, реализуемыми в виде моноатомных
Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 18-07-00778.
© А. Д. Панферов, А. В. Маханьков, 2019
© Саратовский государственный университет им. Н.Г. Чернышевского, 2019 © Программные системы: теория и приложения (дизайн), 2019
или мономолекулярных слоёв на подложках или в свободном состоянии. Графен является одним из таких перспективных материалов [1]. Благодаря чрезвычайно высокой подвижности носителей заряда [2,3] и необычным оптоэлектрическим свойствам [4-6] он может найти применение при разработке новых сверхбыстрых оптоэлектронных устройств [7-9]. Особенно интересны его нелинейные свойства в инфракрасной области спектра и оптическом диапазоне, приводящие к генерации высокочастотных гармоник [10-13]. Данное направление требует сложных экспериментальных технологий и использования уникального оборудования. Теоретические модели происходящих при этом процессов важны как для понимания их физической природы, так и для оценки достижимых параметров и характеристик.
Даже при простоте и прозрачности физических принципов модели, процедура предсказания на её основе наблюдаемых величин может быть вычислительно сложной ресурсоёмкой задачей. Такая проблема рассматривается в настоящей работе, посвященной демонстрации на конкретном примере возможностей разрабатываемого программного комплекса для моделирования поведения носителей заряда в графене в присутствии внешних электрических полей различной природы на основе нового кинетического квантово-полевого подхода [14]. Особенности подхода позволяют рассматривать электрические поля с произвольной зависимостью от времени, в том числе с переменным направлением действия. Единственным ограничивающим условием является требование пространственной однородности поля. Для апробации системы на задаче с реалистическими параметрами выбраны условия экспериментов, результаты которых представлялись в работах [15,16].
1. Физическая модель
Базовое кинетическое уравнение (КУ) для описания электрон-дырочных возбуждений в электрическом поле для Б = 2 + 1 модели графена имеет вид [14]:
(1)
г
¡(Р!,Р2,г) = 1 л (рир2,г) I¿г'А(рьр2,г') [1 - 2/(рьр2,г')]ео8в(г,г'),
г0
где
(2)
Ъ
а параметр
(3)
А (Р1,Р2^}
е«р [Е1Р2 — Е2р1] е2(Р1,Р2,^}
определяется через компоненты двумерного (к = 1, 2} квазиимпульса Рк = Рк — (е/с}Ак (¿} и энергию квазичастиц
(4)
£(#*}= «р ^Р2 = «ру7 (Р1}2 + (Р2}2,
учитывающую специфику закона дисперсии в окрестностях точек Дирака. В этих определениях: Н —приведенная постоянная Планка, е — элементарный заряд, с —скорость света, Vр — скорость Ферми, Рк — компоненты импульса квазичастиц в отсутствии электрического поля, Ек — напряженность и Ак — векторный потенциал электрического поля.
Для численного исследования поведения решений КУ (1) его удобно преобразовать в форму системы из трёх обыкновенных дифференциальных уравнений [17]
(5)
/ = — А и, -1 2
2е 2е
и = А (1—2/} ——V, V = — и,
V Н *
Н
с соответствующими начальными условиями /(¿о} = и(£о} = «(¿о} = 0 для физической системы, находящейся в вакуумном состоянии до включения поля. Ключевым объектом для описания характеристик системы является функция распределения /(Р1,Р2,^}- Функции и(Р1,Р2,г} и «(р1,Р2, ¿} носят вспомогательный характер и описывают особенности эволюции поляризационных эффектов и энергии системы.
Для описания электрического поля коротких оптических импульсов определим его в виде:
(6)
Е1(г} = Е0 С08(^}ехр(- г2/2т2} Е2(г} = Е0 сов(^ + ф} ехр(—¿2/2т2}.
Здесь ш — циклическая частота излучения. Характеристика длительности импульса т может определяться через вспомогательный параметр
а (обычно интерпретируемый как число волн основной частоты в импульсе) соотношением т = а/ш. Сдвиг фаз ф между компонентами вектора поля определяет его поляризацию. В случае ф = 0 поляризация линейна, ф = 0.5п соответствует круговой поляризации.
Значения параметров выбраны в соответствии с условиями эксперимента [16]: ш = 2п • 5,0х1014 Hz, что соответствует длине волны Л = 600nm (оптический диапазон). Амплитуда электрического поля Eo = 2,528х104 V/m. Выбранное значение амплитуды поля соответствует максимальной мощности импульса в 50 pW при фокусировке в пятно диаметром 1,5 pm. Параметр ф не фиксирован. Меняя его, можно моделировать процесс изменения поляризации лазерного импульса.
2. Программная реализация
Реализация процесса моделирования должна обеспечивать получение полной информации о функции распределения носителей на момент прекращения действия электрического поля импульса. Строго говоря это f (pi,p2,t ^ то), а интегрирование необходимо выполнить для —то < t < то. Но, поскольку поле (6) очень быстро убирается экспонентой exp(— t2/2т2), реальный временной интервал t;n <t< tout для определения эволюции функции распределения через систему (5) конечен и составляет несколько т. Его конкретное значение выбирается по результатам численных экспериментов из условия устойчивости результатов интегрирования.
Поскольку в реальных экспериментах измеряются интегральные характеристики функции распределения, необходимо получить достаточно подробную информацию о f(pi,p2,t = tout) = fout(pi,P2) для адекватного воспроизведения значений этих характеристик на последующих этапах моделирования поведения системы.
В общем случае априори ничего не известно о её поведении в импульсном пространстве. По условиям задачи рассматривается только первая зона Бриллюэна, но локализована ли fout(p1,p2) в какой-то относительно небольшой области этой зоны, распределена по нескольким компактным областям или относительно равномерно покрывает всю зону можно узнать только по результатам численных экспериментов. Из этого следует, что обязательным этапом решения задачи будет определение диапазона значений аргументов p1 и p2.
Следующая проблема—выбор шага дискретизации. Если /out (pi,p2) окажется гладкой, можно ограничиться достаточно разреженной сеткой. При сложном характере зависимости функции распределения от pi и Р2 могут потребоваться сетки с очень большим числом узлов. Ниже эти проблемы будут продемонстрированы на конкретных примерах.
Важными упрощающими обстоятельством в рассматриваемом случае является отказ от учета обратного влияния генерируемых носителей на действующее поле и учета диссипативных процессов. Это позволяет решать систему уравнений (5) для каждой точки импульсного пространства независимо и работать с узлами сетки в параллельном режиме. Программно это реализуется с использованием процедур MPI. Выбор MPI определен необходимостью иметь возможность проводить моделирование на массово параллельных системах без ограничения уровня масштабирования.
Для собственно процедуры решения системы ОДУ использована библиотека GSL [18]. Выбор GSL определен наличием в ней объектов различного уровня для организации процедуры решения задачи Коши системы ОДУ и поддержкой разнообразных методов численного решения такой задачи. Имеются реализации нескольких модификаций явного и трех модификаций неявного методов Рунге-Кутта, специальные возможности предусмотрены для решения жестких систем. Эффективность использования различных методов применительно к рассматриваемой задаче была подробно изучена [19].
3. Результаты численных экспериментов
Первая серия рисунков показывает этапы определения вида функции распределения остаточных носителей, образовавшихся в результате действия импульса поля вида (6). Основные его параметры (частота и амплитуда) были определены в разделе 1. Для двух дополнительных параметров выбраны значения: ф = 0 и а = 4. Как было отмечено выше, выбор ф = 0 определяет синфазное поведение x и y компонент поля (6). В этом случае их сумма является линейно поляризованной волной с плоскостью поляризации на диагонали x = y. Выбор а = 4 минимизирует длительность действия поля и затраты времени на выполнение моделирования. Рисунок 1а представляет быстрый обзор зоны Бриллюэна. Он построен по 169 точкам (сетка 13x13) и даёт только приближенное представление об области локализации решения.
5х10-10 4,5х10-10 4х10-10 3,5х10-10 3х10-10 [Г
2х10-10 1,5х10-10 1х10-10 5х10-11 - ¡¡¡¡¡¡Ш Ш ]ШВШ
03 2 чЩЩЁР!!
Р2 0 чЦЩР -2^1 -3^ щдш
Р1
(а) быстрый обзор зоны Бриллюэна
Рисунок 1. Три этапа определения вида функции распределения остаточных носителей
Рисунок 1б построен с учетом полученной информации об области локализации функции распределения в импульсном пространстве. Для него вычисления проводились в 484 точках (сетка 22x22). С учетом уточнения области локализации разрешение по импульсному пространству улучшено в 5 раз. Но рисунок даёт только очень грубое представление о форме распределения и не может гарантировать правильное воспроизведение интегральных характеристик.
На следующей итерации шаг сетки уменьшен ещё в 5 раз по каждой из осей. Результат представлен на рисунке 1в .В этом случае система уравнений (5) решается для 10404 различных наборов параметров (сетка 102x102), что обеспечивает достаточно высокий потенциал для распараллеливания. Полученные данные позволяют быть уверенными в гладкости исследуемого решения и точном воспроизведении его параметров. Кроме того, они наглядно демонстрируют отражение пространственной симметрии действующего электрического поля в импульсном пространстве.
Выбор более реалистичного значения для параметра а =16 уже существенно усложняет задачу. Прежде всего пропорционально увеличивается время действия поля. Но есть и другие факторы, влияющие на вычислительную сложность задачи.
На рисунке 2 а представлены результаты расчётов на сетке, аналогичной использовавшейся для рисунка 1в.
Можно сделать вывод, что функция распределения изменяется.
Примерно на порядок увеличивается её максимальное значение, но при этом «толщина» характерной области локализации существенно уменьшается. В результате сетка 102x102 уже не обеспечивает точного воспроизведения её формы. Для решения этой проблемы число узлов сетки по каждой оси импульсного пространства необходимо увеличить.
На рисунке 2б представлен результат для области 0 < р1 < 0,5, -0,5 < р2 < 0 с шагом сетки, уменьшенным в 4 раза. Можно сделать вывод, что если не предпринимать дополнительных упрощающих действий, то сложность задачи оказывается примерно пропорциональна кубу параметра а. Это может быть существенным фактором, поскольку в реалистических ситуациях этот параметр имеет достаточно большие значения.
(б) сетка 102x102 для области 0 < рх < 0,5, -0,5 < р2 < 0
Рисунок 2. Построение функции распределения для импульса поля со значением параметра а = 16
В рассматриваемой модели есть ещё один свободный параметр манипуляции с которым позволяют варьировать поляризацию падающего излучения. Поведение модели при изменении степени поляризации демонстрируют рисунки 3а и 3б.
(а) ф = 0, 3п
(б) ф = 0, 5п
Рисунок 3. Влияние фазового сдвига ф на форму функции распределения
Параметры поля (за исключением ф) совпадают с использовавшимися при построении первой серии рисунков (рисунок 1). Для построения использовалась сетка значений 102x102. Если на рисунке 1б представлен результат для ф = 0, то на рисунках 3а и 3б показаны
результаты для значений ф = 0,3п и ф = 0,5п соответственно.
Эта серия рисунков показывает изменение симметрии решения. В случае чистой круговой поляризации ф = 0, 5п искомая функция распределения приобретает вращательную симметрию относительно точки pi = p2 = 0. Потенциально это открывает возможность для редуцирования размерности задачи с 2D до 1D с соответствующим уменьшением её вычислительной сложности. Но в общем случае при произвольных значениях ф присутствует только отражательная симметрия относительно прямой pi = p2. При этом само изменение значения ф заметного влияния на время счета не оказывает.
Все представленные выше результаты получены с использованием при решении системы уравнений (5) явного метода, именуемого в документации GSL как «embedded Runge-Kutta Prince-Dormand (8, 9)». Его превосходство было продемонстрировано в работе [19]. Поскольку рассматриваемые варианты реалистических оптических импульсов (6) при больших значениях параметра а сильно отличаются от модели поля, рассматривавшейся в работе [19], было выполнено контрольное сравнение времени решения задачи в расчёте на один узел сетки для нескольких методов, показывавших лучшие результаты. Тесты проводились для значения а = 256 как в области «гребня» функции распределения, так и вдали от него в области фоновых значений, на вычислительных узлах с процессорами Intel® Xeon® E5405 с тактовой частотой 2,0 GHz. В таблице 1 приведены времена счёта для точек в области максимальных значений:
Таблица 1. Время решения задачи для одного узла сетки
Метод интегрирования (по документации GSL) Время (ms)
Explicit embedded Runge-Kutta-Fehlberg (4, 5) Explicit embedded Runge-Kutta Cash-Karp (4, 5) Explicit embedded Runge-Kutta Prince-Dormand (8, 9) A variable-coefficient linear multistep Adams method in Nordsieck form
304,7 194,0 64,3
144,5
Вне этой области времена счёта не превышают приведённых значений.
Из приведённых данных следует, что «embedded Runge-Kutta Prince-Dormand (8, 9)» метод по прежнему справляется с задачей более чем в два раза быстрее ближайшего конкурента.
Заключение
В работе представлены результаты апробации системы моделирования взаимодействия графена с внешними электрическими полями на примере реалистической физической задачи. Рассматривался результат действия на графен короткого оптического импульса. Продемонстрирована возможность детально воспроизводить результат такого воздействия при наличии достаточных вычислительных ресурсов. Особенностью разрабатываемой системы является новая физическая модель, построенная на основе квантового кинетического уравнения [14]. Это позволяет проводить моделирование без ограничений по характеру зависимости действующего поля от времени (его частотных параметров, амплитуды и направления).
Авторы выражают благодарность за поддержку и плодотворные дискуссии С. А. Левенцу, С. А. Смолянскому и С. О. Пирогову.
Список литературы
[1] M. Baudisch, A. Marini, J. D. Cox, T. Zhu, F. Silva, S. Teichmann, M. Massicotte, F. Koppens, L. S. Levitov, F. J. Garcia de Abajo, J. Biegert. "Ultrafast nonlinear optical response of Dirac fermions in graphene", Nature
Communications, 9 (2018), 1018. t34
[2] L. Wang, I. Meric, P. Y. Huang, Q. Gao, Y. Gao, H. Tran, T. Taniguchi, K. Watanabe, L. M. Campos, D. A. Muller, J. Guo, P. Kim, J. Hone, K.L. Shepard, C. R. Dean. "One-dimensional electrical contact to a two-dimensional material", Science, 342 (2013), pp. 614-617.
[3] L. Banszerus, M. Schmitz, S. Engels, M. Goldsche, K. Watanabe, T. Taniguchi, B. Beschoten, C. Stampfer. "Ballistic transport exceeding 28 ^m in CVD
grown graphene", Nano Lett., 16 (2016), pp. 1387-1391. I ' 34
[4] J. Wang, Y. Hernandez, M. Lotya, J. N. Coleman, W. J. Blau. "Broadband nonlinear optical response of graphene dispersions", Advanced Materials, 21 (2009), pp. 2430-2435. 34
[5] Z. Fei, G.O. Andreev, W. Bao, L. M. Zhang, A. S. McLeod, C. Wang, M.K. Stewart, Z. Zhao, G. Dominguez, M. Thiemens, M. M. Fogler, M.J. Tauber, A. H. Castro-Neto, C.N. Lau, F. Keilmann, D.N. Basov. "Infrared nanoscopy of Dirac plasmons at the graphene-SiO2 interface", Nano Lett., 11 (2011), pp. 1701-1705. 34
[6] A. E. Nikolaenko, N. Papasimakis, E. Atmatzakis, Z. Luo, Z. X. Shen, F. Angelis, S.A. Boden, E. Fabrizio, N.I. Zheludev. "Nonlinear graphene metamaterial", Appl. Phys. Lett., 100 (2012), 181109. I ' 34
[7] M. Garg, M. Zhan, T. T. Luu, H. Lakhotia, T. Klostermann, A. Guggenmos,
E. Goulielmakis. "Multi-petahertz electronic metrology", Nature, 538 (2016), pp. 359-363. 34
[8] A. Sommer, E.M. Bothschafter, S.A. Sato, C. Jakubeit, T. Latka, O. Razskazovskaya, H. Fattahi, M. Jobst, W. Schweinberger, V. Shirvanyan, V. S. Yakovlev, R. Kienberger, K. Yabana, N. Karpowicz, M. Schultze,
F. Krausz. "Attosecond nonlinear polarization and light-matter energy transfer in solids", Nature, 534 (2016), pp. 86-90. I 134
[9] T. Higuchi, C. Heide, K. Ullmann, H. B. Weber, P. Hommelhoff. "Light-field-driven currents in graphene", Nature, 550 (2017), pp. 224-228.
t34
[10] N. Kumar, J. Kumar, C. Gerstenkorn, R. Wang, H. Chiu, A. L. Smirl, H. Zhao. "Third harmonic generation in graphene and few-layer graphite films", Phys. Rev. B, 87 (2013), 121406. 134
[11] S. Hong, J. I. Dadap, N. Petrone, P. Yeh, J. Hone, R. M. Osgood, Phys. Rev. X, 3 (2013), 021014. I 34
[12] G.X. Ni, L. Wang, M. D. Goldflam, M. Wagner, Z. Fei, A.S. McLeod, M. K. Liu, F. Keilmann, B. Ozyilmaz, A. H. Castro Neto, J. Hone, M. M. Fogler, D. N. Basov. "Ultrafast optical switching of infrared plasmon polaritons in high-mobility graphene", Nature Photonics, 10 (2016),
pp. 244-247. 34
[13] 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), 041408. I 134
[14] 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", EP.J Web Conf., 138 (2017), 06004.
d t
34 43
[15] I. Gierz, J.C. Petersen, M. Mitrano, C. Cacho, I.C. Edmond Turcu, E. Springate, A. Stohr, A. Kohler, U. Starke, A. Cavalleri. "Snapshots of non-equilibrium Dirac carrier distributions in graphene", Nature Materials, 12 (2013), pp. 1119-1124. 34
[16] K. J. Tielrooij, L. Piatkowski, M. Massicotte, A. Woessner, Q. Ma, Y. Lee, K. S. Myhro, C. N. Lau, P. Jarillo-Herrero, N. F. van Hulst, F. H. L. Koppens. "Generation of photovoltage in graphene on a femtosecond timescale through efficient carrier heating", Nature Nanotechnology, 10 (2015), pp. 437-443.
^34,36
[17] D.B. Blaschke, A. V. Prozorkevich, G. Röpke, C.D. Roberts, S. M. Schmidt, D.S. Shkirmanov, S.A. Smolyansky. "Dynamical Schwinger effect and high-intensity lasers. Realising nonperturbative QED", Eur. Phys. J. D, 55 (2009), 341. d -36
[18] M. Galassi, J. Theiler et all. GNU Scientific Library—GSL 2.5 documentation, https://www.gnu.org/software/gsl/, 1996-2018. 37
[19] С. А. Левенец, Т. Т. Верёвин, А. В. Маханьков, А. Д. Панферов, С. О. Пирогов. «Моделирование динамики безмассовых носителей заряда в двумерной системе», Компьютерные науки и информационные технологии, Материалы международной научной конференции (2 июля-3 июля 2018 года, Саратов, Россия), Издат. центр «Наука», 2018, с. 242-245. © 37 42
Поступила в редакцию 19.12.2018 Переработана 12.02.2019
Опубликована 26.03.2019
Рекомендовал к публикации
к.т.н. С. А. Амелькин
Пример ссылки на эту публикацию:
А. Д. Панферов, А. В. Маханьков. «Моделирование действия коротких оптических импульсов на графен». Программные системы: теория и приложения, 2019, 10:1(40), с. 33-46.
10.25209/2079-3316-2019-10-1-33-46 url http: //pst а. psiras. ru/read/psta2019_l_33-46 .pdf
Эта же статья по-aнглийски:
10.25209/2079-3316-2019-10-1-47-58
Об авторах:
Анатолий Дмитриевич Панферов К.ф.-м.н., зам. начальника ПРЦНИТ Саратовского государственного университета им. Н.Г. Чернышевского по научно-производственной деятельности. Научные интересы: высокопроизводительные вычисления, параллельное программирование, численное решение квантовых кинетических уравнений, моделирование процессов вакуумного рождение частиц в КЭД, генерации носителей в полупроводниках в том числе бесщелевых, процессов на ранних стадиях столкновения релятивистских ядер.
ПИ 0000-0003-2332-0982
e-mail: panferovad@info.sgu.ru
Алексей Владимирович Маханьков
Аспирант Саратовского государственного университета им. Н.Г. Чернышевского. Научные интересы: моделирование физических процессов на высокопроизводительных вычислительных системах, параллельное программирование.
¡з>
e-mail:
0000-0002-9848-9734
next642009@yandex.ru
Sample citation of this publication:
Anatolii Panferov, Alexey Makhankov. "Simulation of the effect of short optical pulses on graphene". Program Systems: Theory and Applications, 2019, 10:1(40), pp. 33-46. (In Russian).
10.25209/2079-3316-2019-10-1-33-46 url; http : //psta. psiras. ru/read/psta2019_l_33-46 .pdf
The same article in English: 10.25209/2079-3316-2019-10-1-47-58
