Hydrodynamic finite-difference time-domain simulation of spatial dispersion and surface modes of thin metal films

Full Text

Введение

Прогресс в производственных и цифровых технологиях и дизайне новых материалов в настоящее время трудно представить без изучения закономерностей и управления взаимодействием фотонов с индивидуальными молекулами среды или различными возбуждениями в твердом теле. Материалы, в которых возбуждаются плазмонные и другие оптические резонансы, позволяют локализовать энергию электромагнитного поля и усилить оптический отклик вещества на несколько порядков по величине. Одними из такими материалов, являются среды с близкой к нулю диэлектрической проницаемостью (НДП, калька с англ. epsilon-near-zero, ENZ), у которых на некоторой частоте (называемой НДП-частотой ${\omega }_{\text{ENZ}}$ ) действительная часть диэлектрической проницаемости равна нулю. Они обладают такими уникальными свойствами, как локализация в своем объеме электрического поля большой напряженности, малая инерция оптического отклика, большая фазовая скорость и малая групповая скорость [1, 2], что может найти применение в интегрированной фотонике [3—6], нелинейной оптике [7—11], оптической микроскопии [1, 2, 11] и даже квантовых технологиях [12]. В этих средах весьма примечательны нелинейно-оптические эффекты генерации высших гармоник (вплоть до девятой) [13] и сверхконтинуума [14], а также гигантское изменение показателя преломления [15] и эффекты преломления света во времени при сверхбыстрых переключениях диэлектрической проницаемости [16]. В отличие от традиционных плазмонных материалов (например золота), материалы с НДП в рабочем диапазоне длин волн обладают гораздо меньшими оптическими потерями [1, 2, 17], что позволяет их рассматривать и как перспективную платформу для реализации нелинейно-оптических преобразований. Более того, существуют материалы, НДП-частоту которых можно настраивать при синтезе [17]. В то же время развитие науки в данных направлениях привело к тому, что моделирование процессов взаимодействия света с веществом стало неотъемлемым шагом в исследовании. В частности, моделирование показывает в таких средах усиление комбинационного рассеяния света и возможность оптического охлаждения твердого тела [18—21].

При моделировании генерации гармоник и других нелинейно-оптических эффектов в средах с НДП развитие получил гидродинамический подход, при котором свободные электроны рассматриваются как заряженная жидкость. Эта модель учитывает конвекцию электронов, гидростатическое давление, магнитные силы. Также эта модель позволяет описывать существенную пространственную дисперсию материала с НДП [22—25], которая связана с гидростатическим давлением заряженной жидкости и сильно проявляется в тонких пленках. Пространственная дисперсия равносильна нелокальности отклика диэлектрической проницаемости, при которой электрическое поле в некоторой области среды приводит к перераспределению заряда в другой части среды. Эта нелокальность приводит к смещению резонансной частоты генерации второй или третьей гармоники в фиолетовую сторону спектра [24], индуцированной анизотропии диэлектрической проницаемости [22], а также к ограничению максимального значения усиления электрического поля внутри материала [26]. Дополнительно к этому, метаматериалы с нулевой эффективной проницаемостью могут обладать наведенной оптической нелокальностью, связанной с вытекающими модами в ограниченных элементах наноструктуры. Связанные электроны и колебательные возбуждения могут быть смоделированы как осцилляторы Лоренца.

Целью данной работы является численное моделирование линейных эффектов пространственной дисперсии сред с НДП с помощью метода конечных разностей во временной области (англ. finite-difference time-domain, FDTD) [27]. На сегодняшний день ни одно существующее коммерческое или свободно распространяемое ПО для моделирования методом FDTD не содержит встроенной модели пространственной дисперсии. Более того, в средах с НДП очень актуально моделирование переходных процессов на временных интерфейсах, для которого необходимы расчеты во временной области. Разработка такой модели позволит в будущем включить в нее и моделировать различные нелинейно-оптические взаимодействия вблизи НДП при распространении коротких лазерных импульсов или при сверхбыстром переключении диэлектрической проницаемости, что является в конечном итоге важным шагом для развития элементной базы фотонных интегральных схем и оптических аналоговых вычислений [28].

Описание метода и параметров моделирования

В работе [29] представлены основные уравнения гидродинамической модели пространственной дисперсии. Запишем эти уравнения в терминах объемной плотности тока проводимости:

$\begin{array}{c}\frac{\partial \overrightarrow{D}}{\partial t}=\text{rot} \overrightarrow{H}-\overrightarrow{J},\\ \rho =\text{div} \overrightarrow{D},\\ \frac{\partial \overrightarrow{J}}{\partial t}+\gamma \overrightarrow{J}={\omega }_{\text{p}}^2{\epsilon }_0\overrightarrow{E}-{\beta }^2\nabla \rho .\end{array}$ (1)

Здесь $\rho$  обозначает индуцированную плотность заряда, ${\omega }_{\text{p}}$  — плазменная частота, $\gamma$  — параметр частоты столкновений электронов с решеткой из модели Друде, учитывающий затухание в среде, ${\beta }^2=3\text{/}5v_{\text{F}}^2$  — коэффициент, появляющийся в выражении для давления идеального ферми-газа. Уравнения (1) записаны в таком виде и порядке для удобства их дискретизации и встраивания в разностную схему FDTD. Как известно [30], в разностной схеме FDTD каждый шаг во времени происходит поочередное обновление электрического и магнитного полей в каждой точке рабочей области. При этом значения электрического поля рассчитываются следующим образом: (i) вначале, интегрируется во времени уравнение Максвелла для ротора магнитного поля, обновляется значение вектора $\overrightarrow{D}$  электрического смещения, (ii) затем находится вектор $\overrightarrow{E}$ , применяя разностную схему для материальных уравнений.

В нашей работе использовался подход, в котором дисперсионный (т. е. зависящий от частоты) отклик среды полностью выражается через вектор $\overrightarrow{J }$ плотности полного тока проводимости [31]. Так как производная плотности тока определяется значением самого электрического поля $\overrightarrow{E}$  (а не производной) и, наоборот, производная вектора $\overrightarrow{D}$  по времени определяется значением плотности тока — то в разностной схеме FDTD значения тока отстоят на полшага по времени от значений $\overrightarrow{E}$  подобно значениям магнитного поля $\overrightarrow{H}$ . Таким образом, во втором порядке точности аппроксимации производных разностная схема приобретает вид (опуская дискретизацию полей в пространстве):

$\begin{array}{c}{\overrightarrow{D}}^{\text{n}+1}={\overrightarrow{D}}^{\text{n}}+\Delta t\cdot \left({\text{rot}}_{\Delta }{\overrightarrow{H}}^{\text{n}+1/2}-{\overrightarrow{J}}^{\text{n}+1/2}\right),\\ {\rho }^{\text{n}+1}={\text{div}}_{\Delta }{\overrightarrow{D}}^{\text{n}+1},\\ \begin{array}{c}{\overrightarrow{J}}^{\text{n}+1/2}=\frac{2-\gamma \Delta t}{2+\gamma \Delta t}{\overrightarrow{J}}^{\text{n}-1/2}+\\ +\frac{2\Delta t}{2+\gamma \Delta t}\left({\omega }_{\text{p}}^2{\epsilon }_0{\text{E}}^{\text{n}}-{\beta }^2{\text{grad}}_{\Delta }{\rho }^{\text{n}}\right),\end{array}\end{array}$ (2)

где $\Delta t$  — шаг по времени, верхними индексами указаны номера шагов, символами ${\text{rot}}_{\Delta }$ , ${\text{div}}_{\Delta }$ , ${\text{grad}}_{\Delta }$  обозначены пространственные конечно-разностные аппроксимации дифференциальных операторов ротора, дивергенции и градиента соответственно. Вектор $\overrightarrow{D}$  в данной схеме содержит вклад только от связанных электронов. Так как данная модификация разностной схемы сопряжена с введением дополнительных пространственных производных (а именно градиента плотности заряда), то реализовать ее в имеющихся пакетах для моделирования нет возможности. Поэтому алгоритм FDTD в двумерной ТМ-формулировке (с векторами E_x, E_z, H_y) с применением данной разностной схемы для пространственной дисперсии был реализован самостоятельно.

Геометрия моделируемой системы показана на рис. 1. Тонкая металлическая пленка толщиной d = 20 нм, помещена на стеклянную подложку с показателем преломления 1.5 так, что границы пленки совпадают с плоскостью XY. Рабочая область моделирования составляла 5200×400 нм² с размером ячейки равномерной вычислительной сетки 4 нм. Граничными условиями выбраны идеально сочетающиеся слои (PML) в количестве 32 слоев в направлении оси X и 16 слоев в направлении оси Z со стандартными профилем в формулировке растяжения координат (stretched-coordinate). Источник представлял собой точечный магнитный диполь, расположенный в середине металлической пленки, и направление которого перпендикулярно плоскости моделирования. Импульс излучения диполя имел профиль огибающей во времени в виде первой производной функции Гаусса шириной 10 шагов во времени и центральную частоту несущей 500 ТГц. Моделирование осуществлялось 2000 шагов.

 

Рис. 1. Схематическое изображение моделируемой системы, состоящей из тонкой металлической пленки толщиной d, расположенной на подложке из SiO₂ и окруженной воздухом. Плоскость рисунка совпадает с плоскостью XZ. Распространяющий плазмон-поляритон имеет ТМ-поляризацию (вектор E лежит в плоскости рисунка). Внутри пленки с пространственной дисперсией электрическое поле не перпендикулярно волновому вектору.

 

Параметры материала: проницаемость связанных электронов ε_b = 5.4, плазменная частота ω_p = 1.38·10¹⁶ рад/с, γ = 25·10¹² рад/с, скорость Ферми v_F = 3/137скорости света. В ходе моделирования сохранялись значения поля внутри металлической пленки [1].

Результаты и их обсуждение

Проводилось моделирование собственных мод тонкой металлической пленки в двух случаях — отсутствия и наличия у пленки пространственной дисперсии. Для построения дисперсионных зависимостей осуществлялось двумерное пространственно-временное преобразование Фурье. На рис. 2 представлены результаты моделирования в виде дисперсионных диаграмм. Хорошо заметны поверхностные моды, локализованные внутри и вблизи тонкого металлического слоя. Как известно, в таких пленках вблизи НДП-частоты могут возбуждаться поверхностные моды, которые называют модами Берремана для области волновых чисел внутри светового конуса, и модами вблизи НДП, если волновой вектор лежит вне светового конуса [32]. Наклон дисперсионных кривых этих мод и величина нормальной компоненты электрического поля зависят от толщины пленки, и при достаточно малой величине последней наблюдается усиление нелинейно-оптических эффектов [33]. Оба этих явления, а именно малая групповая скорость и усиление нормальных к поверхности компонент поля, предлагаются в качестве механизмов усиления [34]. Асимптоты дисперсионных кривых высокочастотной и низкочастотной мод ниже плазменной частоты также являются постоянными функциями на бесконечности, причем дисперсионная кривая высокочастотной моды в общем случае не монотонна [32].

 

Рис. 2. Дисперсионные диаграммы собственных мод тонкой металлической пленки в случае отсутствия (а) и наличия (б) пространственной дисперсии, полученные с помощью моделирования. Параметры пленки: d = 20 нм, ε_b = 5.4, ω_p = 1.38·10¹⁶ рад/с, γ = 25·10¹² рад/с, скорость Ферми v_F = 3/137 скорости света. Штриховой линией обозначено положение частоты НДП, штрихпунктирной линией обозначена кривая, определяющая положение моды Берремана (внутри светового конуса) и НДП-моды (вне светового конуса) в случае наличия пространственной дисперсии (б).

 

Тем не менее, в случае наличия у материала пространственной дисперсии (рис. 2б) видно, что дисперсионные кривые поверхностных мод монотонно возрастают. Более того, их асимптоты имеют положительный наклон. Это означает, что групповая скорость НДП-моды перестает быть малой, что ограничивает усиление нелинейно-оптического отклика. Кроме того, как и другие моды, связанные с какой-либо неоднородностью, их невозможно возбудить падающим светом непосредственно — для возбуждения НДП-моды необходимо выполнение специальных условий фазового синхронизма или согласование поверхностных импедансов.

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

Для проверки полученных результатов, дисперсионные соотношения были аналитически найдены из уравнений Максвелла. Дисперсионное уравнение любой моды имеет вид g(ω, k_x) = 0, т. е. дисперсионные кривые соответствуют нулю функции g(ω, k_x). Поэтому для визуализации дисперсионных кривых был построен натуральный логарифм функции g(ω, k_x)^-1, они соответствуют полюсам этой функции. Аналитический вид функции g(ω, k_x) для металлического слоя, ограниченного с обеих сторон диэлектриками, может быть найден в работах [18, 32]:

$g\left(\omega ,k_{\text{x}}\right)=i\text{tg}\left(k_{2\text{z}}d\right)\left(\frac{{\kappa }_1{\kappa }_3}{{\kappa }_2}+{\kappa }_2\right)-\left({\kappa }_1+{\kappa }_3\right),$ (3)

где ${\kappa }_{\text{j}}=k_{\text{jz}}\text{/}{\epsilon }_{\text{j}}=\sqrt{\left({\epsilon }_{\text{j}}{\omega }^2\text{/}{c}^2-k_{\text{x}}^2\right)}\text{/}{\epsilon }_{\text{j}}$ , индекс $j=\mathrm{1,2,3}$  для воздуха, металла и стекла соответственно, в случае металла ${\epsilon }_2={\epsilon }_{\text{b}}-{\omega }_{\text{p}}^{\text{2}}/\left({\omega }^2+i\omega \gamma \right)$ , d — толщина металлической пленки. На рис. 3а показаны дисперсионные кривые собственных поверхностных мод плазмон-поляритонов для среды без пространственной дисперсии. Как видно в целом поведение дисперсионных кривых соответствует результатам численного моделирования, что говорит о корректности численного метода. Незначительные отклонения для мод выше плазменной частоты объясняются тем, что при моделировании записывалось поле только внутри металлического слоя, тем самым моды, локализованные преимущественно за пределами слоя, оказались менее проявлены.

 

Рис. 3. Дисперсионные диаграммы (а) мод тонкой металлической пленки в случае отсутствия пространственной дисперсии и (б) поверхностного плазмон- поляритона, возбуждаемого на поверхности материала с пространственной дисперсией. Скорость Ферми v_F = 24/137 скорости света. Штриховой линией обозначено положение частоты НДП, которая определяет также положение моды Берремана (внутри светового конуса) и НДП-моды (вне светового конуса) при отсутствии пространственной дисперсии (а).

 

В случае наличия пространственной дисперсии функция g(ω, k_x) была найдена, решая уравнения Максвелла совместно с материальным уравнением (1) гидродинамической модели газа свободных электронов:

$g\left(\omega ,k_{\text{x}}\right)={\kappa }_1-{\kappa }_2+\frac{{\epsilon }_2-{\epsilon }_{\text{b}}}{{\epsilon }_{\text{b}}{\epsilon }_2}\frac{k_{\text{x}}^2}{k_{\text{z}}^{\parallel }}$, (4)

где $k_z^{\parallel }=\sqrt{\left({\epsilon }_2/{\epsilon }_{\text{b}}\left({\omega }^2+i\omega \gamma \right)/{\beta }^2-k_{\text{x}}^2\right) }$ находится из условия ε_‖(ω, k^‖) = 0. При $\beta \to 0$  уравнение (4) переходит в дисперсионное соотношение для обычных поверхностных плазмонов. На рис. 3б построены дисперсионные кривые поверхностного плазмон-поляритона на одиночной границе материала с пространственной дисперсией. Для того чтобы подчеркнуть вклад от пространственной дисперсии, параметр β был искусственно увеличен в 8 раз по сравнению со случаем численного моделирования. Это говорит о том, что эффекты пространственной дисперсии сильнее проявляются в ограниченных материалах. Причина изменения наклона НДП-моды лежит в том, что пространственная дисперсия материала — следствие возбуждения в газе свободных электронов акустических волн. Плазменные колебания благодаря этому перестают быть стоячими. С данными колебаниями плотности заряда связана продольная электрическая волна, распространяющаяся в объеме среды со скоростью  $\beta$. На границах материала происходит смешивание обычной поперечной световой волны и продольной плазменной волны, так как граничные условия выполняются для полного электрического поля. Поэтому чем больше пространственных границ у материала, тем больше происходит гибридизация продольной и поперечных волн. Заметим также, что сдвиг высокочастотной моды Берремана оказывается больше для металлической пленки по сравнению с одиночной границей полубесконечного материала.

Заключение

Таким образом, реализована гидродинамическая модель пространственной дисперсии материала при моделировании методом FDTD. С помощью данного численного метода были найдены дисперсионные кривые собственных поверхностных мод тонкой металлической пленки в случаях с пространственной дисперсией и без нее. Для проверки корректности метода дисперсионное уравнение также было найдено аналитически. В тонких пленках из материала с малым показателем преломления, вблизи частоты НДП существуют мода НДП и мода Берремана. В ходе моделирования было выявлены изменение наклона дисперсионных кривых НДП-моды тонкой пленки и сдвиг частоты моды Берремана в синюю сторону оптического спектра. Первый эффект приводит к увеличению групповой скорости, что негативно отражается на величине нелинейно-оптических эффектов. Это ограничивает величину усиления электрического поля в объеме материала с близкой к НДП. Также пространственная дисперсия приводит к сдвигу частоты моды Берремана, что должно быть учтено при разработке компактных нелинейно-оптических устройств преобразования частоты света, использующих данные материалы. Моделирование показало, что в тонких пленках происходит возрастание величины данных эффектов. Большее пространственное ограничение материала с одной стороны приводит к возрастанию нелинейно-оптического отклика материала, с другой стороны — к более заметному проявлению пространственной дисперсии. Эти особенности должны быть учтены при проектировке соответствующих устройств или подготовке экспериментов. Данное исследование и представленный алгоритм является шагом на пути масштабного распространения материалов с НДП для их использования в нелинейно-оптических устройствах фотонных интегральных схем.

Работа финансировалась за счет средств программы стратегического академического лидерства «Приоритет 2030» Казанского федерального университета при Правительстве Российской Федерации.

About the authors

А. R. Gazizov

Kazan (Volga Region) Federal University

Author for correspondence.
Email: almargazizov@kpfu.ru

Институт физики

Russian Federation, Kazan

E. A. Izbasarova

Kazan (Volga Region) Federal University

Email: almargazizov@kpfu.ru

Институт физики

Russian Federation, Kazan

References

Supplementary files