Известия РАН. Серия физическая, 2020, T. 84, № 3, стр. 418-424

ВВЕДЕНИЕ

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

В настоящей работе рассматриваются вопросы разработки и реализации математических моделей для проектирования плазмонных схем на основе полупроводниковых квантовых точек и наноструктурированного графена. Для этого, была создана база данных характеристик КТ на основе соединений A3B5 и A2B6, вмещающая в себя информацию о параметрах энергий потолка валентной зоны, дна зоны проводимости, ширины запрещенной зоны и т.д. Разработанные математические методы позволили производить автоматизированный расчет энергетической структуры уровней, частот и дипольных моментов, скоростей релаксации и др.

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

МОДЕЛЬ РАСЧЕТА ЗАВИСИМОСТЕЙ ПОЛОЖЕНИЯ ЭНЕРГЕТИЧЕСКИХ УРОВНЕЙ ДЛЯ ПОЛУПРОВОДНИКОВЫХ КТ ЗАДАННОГО РАЗМЕРА

Для расчета размерных зависимостей ограниченного числа энергетических уровней дырки в валентной зоне КТ может быть использовано оценочное соотношение [3–6]

где ${{E}_{V}}$ – энергия для потолка валентной зоны КТ (выраженной в эВ); $e = 1.6 \cdot {{10}^{{ - 19}}}$ – заряд электрона (${\text{Кл}}$); ${{m}_{h}}$ – эффективная масса дырки;$~a~$ – радиус КТ (м);$~\hbar = 1.05 \cdot {{10}^{{ - 34}}}\,\,{\text{Дж}} \cdot {\text{с}}~$ – постоянная Планка; $~n$ – главное квантовое число для дырки; $l$ – орбитальное квантовое число для дырки, ${{\chi }_{{nl}}}$ – корни сферического уравнения Бесселя 1-го рода. Для расчета размерных зависимостей энергетических уровней электрона в зоне проводимости КТ может быть использовано оценочное соотношение [5, 6]: где ${{E}_{c}}$ – энергия дна зоны проводимости КТ (эВ); ${{m}_{e}}$ – эффективная масса электрона; $n{\kern 1pt} '$ – главное квантовое число для электрона, $l{\kern 1pt} '$ – орбитальное квантовое число для электрона; ${{\chi }_{{n{\kern 1pt} 'l{\kern 1pt} '}}}$ – корни сферического уравнения Бесселя 1-го рода. Для расчета частоты внутризонного дырочного перехода может быть использовано соотношение [7]: где $m,~n$ и $k$ задаются произвольным образом в рамках пользовательских диапазонов, параметр $l$ подчиняется правилу отбора $l = k \pm 1.$ На основе частоты перехода может быть рассчитана длина волны межуровневого перехода в виде:$~{{\lambda }_{{mn}}} = \frac{{2\pi c}}{{{{\omega }_{{mn}}}}}.$ Для расчета частоты внутризонного (электронного) перехода может быть использовано соотношение [7]: где $m{\kern 1pt} ',~n{\kern 1pt} '$ и $k{\kern 1pt} '$ задаются произвольным образом в рамках интересующих диапазонов, параметр $l{\kern 1pt} '$ подчиняется правилу отбора $~l{\kern 1pt} ' = k{\kern 1pt} ' \pm 1.$ На основе частоты перехода может быть рассчитана длина волны межуровневого перехода в виде ${{\lambda }_{{m{\kern 1pt} 'n{\kern 1pt} '}}} = \frac{{2\pi c}}{{{{\omega }_{{m{\kern 1pt} 'n{\kern 1pt} '}}}}}.$ Тогда, для расчета частоты межзонного перехода используется соотношение [7]:

Правило отбора в данном случае имеет вид $k{\kern 1pt} ' - l = 0.$ Длина волны соответствующего перехода имеет вид $~~{{\lambda }_{{m{\kern 1pt} 'n{\kern 1pt} }}} = \frac{{2\pi c}}{{{{\omega }_{{m{\kern 1pt} 'n}}}}}.$ Дипольный момент внутризонного перехода в КТ, согласно [8], может быть рассчитан следующим образом:

При переходе к тройному интегралу в сферических координатах данный параметр может быть рассчитан так:

где $\varepsilon $ – диэлектрическая проницаемость КТ, $r$ – радиус-вектор, определяющий расстояние носителя заряда до центра КТ, $\theta $ – азимутальный угол, $~\varphi $ – полярный угол, ${{\psi }_{i}},~$ ${{\psi }_{f}}$ – волновые функции уровней, с которого и на который осуществляется переход соответственно; $~R$ – радиус КТ; $\psi _{i}^{*}\left( {r,{\theta },{\varphi }} \right)$ – сопряженная волновая функция уровня, с которого осуществляется переход. Согласно [9] волновая функция для сферической КТ определяется следующим образом:

либо

где ${{j}_{l}}\left( x \right),$ ${{j}_{{l + 1}}}\left( x \right)$ – сферические функции Бесселя 1-го рода и l-го и l + 1-го порядков соответственно; ${{J}_{{l + 1{\text{/}}2}}}\left( x \right),$ ${{J}_{{l + 3{\text{/}}2}}}\left( x \right)$ – функции Бесселя 1-го рода l + 1/2-го и l + 3/2-го порядков соответственно; ${{Y}_{{lm}}}\left( {{\Theta },\phi } \right)$ – сферические функции; ${{k}_{{nl}}} = \frac{{{{\chi }_{{nl}}}}}{R};$ ${{\chi }_{{nl}}}$ – n-ый корень функции Бесселя l-го порядка; $n,l,m$ – главное, орбитальное, магнитное квантовые числа частицы, соответственно. Подставляя (8) для двух различных уровней в (7), получим: где $n{\kern 1pt} ',$ $l{\kern 1pt} ',$ $m{\kern 1pt} '~~$ и $n,$ $l,~$ $m$ – главное, орбитальное и магнитное квантовые числа уровней, между которыми осуществляется переход. В общем виде сферические функции представляют собой следующие выражения: где $\phi $ – азимутальный угол; $\theta $ – телесный угол; ${{{\Theta }}_{{lm}}}\left( \theta \right)$ – радиальные функции, имеющие вид: где $P_{l}^{m}\left( {\cos \left( \theta \right)} \right)$ – присоединенные полиномы Лежандра (см. [10]) в виде:

При $~m = 0~$ полином $P_{l}^{m}$ совпадает с полиномом Лежандра ${{P}_{l}}$ вида:

Далее мы полагаем, что магнитное квантовое число равно m = 0, тогда:

Дипольный момент межзонного перехода может быть приблизительно оценен исходя из соотношения [11]:

где ${{\omega }_{{m{\kern 1pt} 'n}}}$ – частота перехода; ${{E}_{{m{\kern 1pt} 'n}}}$ – расстояние между уровнями $~m{\kern 1pt} '$ и $n,$ $\Delta {{~}_{0}}$ – спин-орбитальное расщепление. Для расчета энергий, частот и дипольных моментов переходов в КТ, использовались формулы (1)–(14), объединенные в единый программный модуль [12]. В частности, результаты расчетов зависимостей энергии межзонного перехода в диапазоне радиусов от 3 до 6 нм для КТ на основе п/п CdSe представлены на рис. 1. Полученная с помощью аналитических выражений зависимость совпадает с той, что была получена в работе [6]. Однако для КТ с радиусом менее 3 нм будет наблюдаться расхождение с результатами наилучшей аппроксимации из работы [6]. Это связано с необходимостью учета эффектов кулоновского взаимодействия в режиме сильного конфаймента [9] при условии $a < {{a}_{{ex}}},$ где ${{a}_{{ex}}}$ – боровский радиус экситона.

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ПРОВОДИМОСТИ ГРАФЕНА И МЕТОД КОНЕЧНЫХ РАЗНОСТЕЙ ВО ВРЕМЕННОЙ ОБЛАСТИ

В качестве основы для формирования межузловых соединений в плазмонных схемах могут быть выбраны как металлические, так и более перспективные, графеновые волноводы. Графен [13] является двумерным материалом на основе гексагональной решетки атомов углерода толщиной в один атом. Уникальной особенностью такой структуры является возможность наблюдения высокотемпературной проводимости как в чистом материале [14], так и его допированных модификациях [15], что может сделать графен основой интегральных схем будущего. Преимущество графена как основы для плазмонных схем заключается в возможности достижения предельной локализации электромагнитной волны на масштабах в несколько нм [16]. Электропроводность графена на частоте $\omega $ может быть описана на основе формализма Кубо [17] и определяется выражением:

где ${{\mu }_{с}}$ – химический потенциал, $1{\text{/}}\tau $ – скорость рассеяния носителей заряда (электронов), T = = 300 К – температура, k – постоянная Больцмана. Первое слагаемое в уравнении (16) описывает вклад внутризонной проводимости, а второе – межзонной. Формирование ППП на поверхности графен–диэлектрик возможно в случае $\hbar \omega < 2{{\mu }_{с}},$ когда мнимая часть проводимости ${\text{Im}}\left( {\sigma \left( {\omega ,{{\mu }_{с}},\tau ,T} \right)} \right)$ переходит от положительных к отрицательным значениям. В частности, для того чтобы ППП существовали при 1.55 мкм, химический потенциал графена должен составлять не менее 0.5 эВ. В случае представления проводимости в виде $\sigma \left( \omega \right)$ = ${{\sigma }_{r}}\left( \omega \right) + i{{\sigma }_{i}}\left( \omega \right)$ диэлектрическая проницаемость графена задается в виде: где ${{\varepsilon }_{r}}\left( \omega \right)$ – диэлектрическая проницаемость диэлектрика, в котором находится лист графена, ${\Delta }$ – толщина графена (обычно, выбирают ${\Delta }$ ~ ~ 1 нм). Таким образом, формирование ППП в листе графена возможно лишь при условии ${{\varepsilon }_{r}}\left( \omega \right) - {{{{\sigma }_{i}}\left( \omega \right)} \mathord{\left/ {\vphantom {{{{\sigma }_{i}}\left( \omega \right)} {\left( {{{\varepsilon }_{0}}\omega {\Delta }} \right)}}} \right. \kern-0em} {\left( {{{\varepsilon }_{0}}\omega {\Delta }} \right)}} < 0.$ Относительная диэлектрическая проницаемость листа графена в вакууме определяется следующим уравнением [1]: где ${{\sigma }_{0}} = {{{{e}^{2}}kT\tau } \mathord{\left/ {\vphantom {{{{e}^{2}}kT\tau } {\left( {\pi {{\hbar }^{2}}{\Delta }} \right)}}} \right. \kern-0em} {\left( {\pi {{\hbar }^{2}}{\Delta }} \right)}}\left( {{{{{\mu }_{с}}} \mathord{\left/ {\vphantom {{{{\mu }_{с}}} {\left( {kT} \right)}}} \right. \kern-0em} {\left( {kT} \right)}} + 2{\text{ln}}\left( {{{e}^{{{{ - {{\mu }_{с}}} \mathord{\left/ {\vphantom {{ - {{\mu }_{с}}} {\left( {kT} \right)}}} \right. \kern-0em} {\left( {kT} \right)}}}}} + 1} \right)} \right)$ и ${{\varepsilon }_{{0~}}}$ – диэлектрическая проницаемость воздуха. При этом формируемая в одиночном листе ППП-волна характеризуется следующим значением волнового вектора: где ${{\tilde {\sigma }}_{{intra}}} = {{{{\sigma }_{0}}{\Delta }} \mathord{\left/ {\vphantom {{{{\sigma }_{0}}{\Delta }} {\left( {1 + i\omega \tau } \right)}}} \right. \kern-0em} {\left( {1 + i\omega \tau } \right)}}.$ При выборе значений ${{\mu }_{с}} = 0.5\,\,{\text{эВ}},$ $\tau = ~3.85\,\,{\text{пс}}$ и частоте источника 3 · · 10¹² Гц длина волны ППП составит λ_ППП ≈ ≈ 51 мкм, а эффективная длина пробега $~{{L}_{{{\text{ППП}}}}}$ = = ${{{{\lambda }_{0}}} \mathord{\left/ {\vphantom {{{{\lambda }_{0}}} {\left( {4\pi {\text{Im}}\left( {{{{{\beta }_{{{\text{ППП}}}}}} \mathord{\left/ {\vphantom {{{{\beta }_{{{\text{ППП}}}}}} {{{k}_{0}}}}} \right. \kern-0em} {{{k}_{0}}}}} \right)} \right)}}} \right. \kern-0em} {\left( {4\pi {\text{Im}}\left( {{{{{\beta }_{{{\text{ППП}}}}}} \mathord{\left/ {\vphantom {{{{\beta }_{{{\text{ППП}}}}}} {{{k}_{0}}}}} \right. \kern-0em} {{{k}_{0}}}}} \right)} \right)}}$ составит 396 мкм.

Для численного моделирования формирования ППП на графене использовался метод FDTD для уравнений Максвелла вида [1]:

где $\varepsilon _{r}^{*}$ и ${{c}_{0}}$ – относительная диэлектрическая проницаемость среды и скорость света, D – плотность потока, E – напряженность электрического поля, H – напряженность магнитного поля, μ – магнитная проницаемость среды. При переходе к нормированным параметрам $\tilde {E} = \sqrt {{{{{\varepsilon }_{0}}} \mathord{\left/ {\vphantom {{{{\varepsilon }_{0}}} {{{\mu }_{0}}}}} \right. \kern-0em} {{{\mu }_{0}}}}} E$ и $\tilde {D} = {1 \mathord{\left/ {\vphantom {1 {\sqrt {{{\varepsilon }_{0}}{{\mu }_{0}}} D}}} \right. \kern-0em} {\sqrt {{{\varepsilon }_{0}}{{\mu }_{0}}} D}}$ уравнения Максвелла перепишутся в следующей форме:

Уравнения Максвелла для поперечной (ТЕМ) моды, когда поляризованная волна распространяется в направлении z, имеют вид:

Дискретизация уравнения Максвелла с использованием метода центральной производной и шагом ${\Delta }t = {{{\Delta }z} \mathord{\left/ {\vphantom {{{\Delta }z} {2{{c}_{0}}}}} \right. \kern-0em} {2{{c}_{0}}}},$ где ${\Delta }z$ – пространственный размер (рис. 2), дает:

Электрическое поле вычисляется по формуле (21в). Объединив уравнения (18) и (21в), получим:

где $B = {{\tau {{\sigma }_{0}}} \mathord{\left/ {\vphantom {{\tau {{\sigma }_{0}}} {{{\varepsilon }_{0}}}}} \right. \kern-0em} {{{\varepsilon }_{0}}}}$. Уравнение (24) зависит от частоты, которую можно преобразовать во временную зависимость и получить набор уравнений: где, $M\left( k \right)$ = ${{{\Delta }t{{\sigma }_{0}}} \mathord{\left/ {\vphantom {{{\Delta }t{{\sigma }_{0}}} {{{\varepsilon }_{0}}}}} \right. \kern-0em} {{{\varepsilon }_{0}}}} - {{{\Delta }tB} \mathord{\left/ {\vphantom {{{\Delta }tB} \tau }} \right. \kern-0em} \tau },$ $N\left( k \right) = {{{\Delta }t{{\sigma }_{0}}} \mathord{\left/ {\vphantom {{{\Delta }t{{\sigma }_{0}}} {{{\varepsilon }_{0}}}}} \right. \kern-0em} {{{\varepsilon }_{0}}}}$ и $L\left( k \right) = {{{\Delta }tB} \mathord{\left/ {\vphantom {{{\Delta }tB} \tau }} \right. \kern-0em} \tau },$ $i,$ $j$ – индексы по пространственным координатам, $n$ – индекс по временной координате, ${\Delta }x$ – шаг по пространственной оси, ${\Delta }t$ – шаг по временной оси. После приведения дискретные уравнения Максвелла примут вид:

Предварительное моделирование формирования и распространения поверхностных электромагнитных волн проводилось с использованием в качестве источника магнитного диполя в виде:

Изначально была проведена калибровка реализованного метода FDTD по известным результатам работ [1, 2, 16]. Далее было проведено полное моделирование поля с источником (27) на длинах волн 2, 4 и 8 мкм для структуры из двух графеновых листов и показано формирование ППП в системе (рис. 3). Для удобства было разработано онлайн-приложение, в котором на стороне пользователя находится построитель макета двумерных плазмонных схем, в то время как блок численного расчета и моделирования электромагнитного поля реализован на стороне сервера. Разработанные подходы и их реализация в настоящем варианте не позволяют проводить полномасштабное моделирование процессов резонансного взаимодействия КТ и ППП в графене, т.е. моделирование для нескольких КТ вблизи графена, что является темой последующих исследований.

Вместе с тем, версия программы, адаптированная для работы с тонкими пленками произвольных проводящих материалов, позволила провести полное полевое моделирование процессов формирования и выявить особенности управления поверхностными плазмон-поляритонами в планарных устройствах на основе наноструктурированных пленок золота. В перспективе запланировано расширение программы для реализации трехмерного моделирования и анализа как отдельных локализованных плазмон-поляритонных структур [18], так и целых массивов на их основе [19].

ЗАКЛЮЧЕНИЕ

В работе рассмотрены математические модели и численные алгоритмы, необходимые для проектирования функциональных устройств на основе двумерных материалов и полупроводниковых квантовых точек. В частности, разработаны модели для расчета основных оптических характеристик КТ на основе полупроводниковых соединений A3B5 и A2B6. Составлена и программно реализована численная модель метода конечных разностей во временной области для моделирования распространения плазмон-поляритонных мод в плазмонных структурах на основе графена. Проведена проверка корректности работы алгоритмов на основе сравнения результатов симуляции с данными экспериментальных и теоретических работ.

Работа выполнена в рамках договора 2226ГС1/ 37022 с фондом содействия инновациям и поддержана грантом Российского научного фонда № 19-72-20154.