Журнал вычислительной математики и математической физики, 2019, T. 59, № 1, стр. 118-127

1. ВВЕДЕНИЕ

Один из наиболее интересных и актуальных разделов нанофотоники посвящен изучению физических явлений, связанных с возможностью управления световым излучением и получения сильных сфокусированных электромагнитных полей в масштабах нескольких атомов [1]–[3]. Такая фокусировка достигается посредством использования эффекта плазмонного резонанса (ПР). Важнейшим свойством плазмонов является их способность манипулировать светом, а именно – направлять его вдоль металлических полос или фокусировать его в отверстиях, сужениях, зазорах между частицами или вблизи участков поверхностей с резким изменением профиля. В этом случае концентрация поля происходит в объемах, имеющих размеры в сотни раз меньше длины волны излучения, но интенсивность электрического поля может достигать нескольких порядков величин [3], [4]. Практика моделирования поведения ПР основана на классической электромагнитной теории Максвелла, в которых коллективное поведение электронов описывается в рамках эмпирических диэлектрических постоянных [1], [2], а распределение электромагнитного поля описывается решением системы уравнений Максвелла. Подобный подход предсказывает монотонное возрастание электрического поля по мере уменьшения зазора между частицами, что побуждает исследователей к разработке нанотехнологий для практической реализации плазмонных структур с субнаноразмерными зазорами [5], [6].

Кроме прочего, кластеры плазмонных частиц с субнанометровыми зазорами оказываются весьма привлекательными в качестве составляющих элементов плазмонных систем, поскольку они позволяют манипулировать положением ПР и его амплитудой в широком диапазоне длин волн. В ряде работ было показано, что квантовомеханические эффекты, такие как нелокальное экранирование и туннелирование электронов, приобретают существенное значение, по мере того, как размеры зазора приближаются к шкале субнанометра. Учет возникающих при этом квантовых эффектов в значительной степени разрушает классическое описание поведения полей в подобных системах. Вместе с тем недавние теоретические и экспериментальные исследования [7]–[9] показывают, что по мере того, как размер зазора попадает в нанометровую и субнанометровую шкалу, квантовая природа электронов и нелокальное экранирование могут значительно изменять плазмонный отклик рассматриваемой структуры [10], [11]. В этом квантовом режиме оказывается, что классическое описание не учитывает фактическую локализацию поверхностных зарядов, вызванных внешним электромагнитным полем [12], [13]. Сдвиг наведенных поверхностных зарядов относительно геометрических границ металла тесно связан явлением нелокальности и приводит к экранированию зарядов внутри структуры [14], [15]. Этот эффект может быть описан с использованием концепции нелокальных диэлектрических функций с различными уровнями сложности в имеющихся теоретических построениях [16], [17].

Можно схематически проиллюстрировать влияние квантовых эффектов на электромагнитные свойства ПР наноструктур с различными зазорами. Размеры зазора до 2–5 нм (в зависимости от материала структуры) соответствуют классическому режиму, для которого система уравнений Максвелла правильно описывает смещение ПР в длинноволновую область и величину его амплитуды при уменьшении зазора. Как только размер зазора становится меньше 2–5 нм, система переходит в квантовый режим, что требует более детального исследования. В данном диапазоне эффект нелокального экранирования (ЭНЭ) выступает как доминирующий квантовый эффект. В этом случае поведение ПР качественно сходно с предсказаниями локальной классической модели, но возникают существенные количественные различия как в амплитуде ПР, так и его положении. Таким образом, для правильного количественного описания происходящих процессов и связанных с ними изменений в плазмонном отклике в диапазоне размеров 0.5–2 нм необходимо существенно учитывать ЭНЭ [7], [17]. По мере дальнейшего уменьшения зазоров <0.5 нм преобладающий вклад вносит эффект туннелирования электронов, который полностью изменяет поведение плазмонного отклика [18]–[20].

В данной работе метод Дискретных источников (МДИ) модифицируется применительно к исследованию линейного кластера плазмонных частиц с наноразмерным зазором, обеспечивающим переход в квантовый режим, с учетом ЭНЭ [7], [17]. Одним из проявлений ЭНЭ является появление продольных электромагнитных полей. Когда размер элементов структуры становится существенно меньше свободного пробега электронов в веществе, внутри формируется объемный заряд, который под действием внешнего поля образует ток. В этом случае внутреннее электрическое поле Е перестает быть чисто поперечным (div E = 0), и для адекватного описания происходящих процессов возникает необходимость привлечения дополнительно продольных полей (rot E = 0) [21]. При этом учет наличия продольных полей может осуществляться на основе различных теорий, в данной работе учет осуществляется на основе наиболее востребованной в настоящее время модели обобщенного нелокального отклика (ОНО) (Generalized Nonlocal Optical Response (GNOR)) [17], [22], [23]. В рамках GNOR проводится обобщение гидродинамической теории нелокального отклика посредством учета диффузии электронов [23]. Именно это обстоятельство позволило модели GNOR оказаться наиболее подходящим инструментом описания процессов, происходящих в кластере наночастиц с квантовым зазором [24]. На основе компьютерной модели МДИ проведен сравнительный анализ классической и модели ОНО, по их влиянию, как на величину ПР, так и его положение в спектральной области. Исследовано поведение интенсивности поля в зазоре между частицами в спектральном диапазоне.

2. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ОБОБЩЕННОГО НЕЛОКАЛЬНОГО ОТКЛИКА

Сразу отметим, что введение продольных компонент поля для описания ЭНЭ в рамках ОНО проводится на основе обобщения закона Ома для тока проводимости внутри [23], то есть с осуществлением перехода следующего вида:

В результате чего меняется соответствующее уравнение системы Максвелла для магнитного поля.

Будем рассматривать рассеяние поля электромагнитной плоской волны $\{ {{{\mathbf{E}}}_{0}},\;{{{\mathbf{H}}}_{0}}\} $ линейным кластером, состоящим из двух проницаемых осесимметричных частиц с внутренними областями${{D}_{{1,2}}}$, расположенных в ${\text{ }}{{R}^{3}}$. Пусть частицы имеют гладкие поверхности $\partial {{D}_{{1,2}}} \in {{С }^{{(2,\alpha )}}}$ (пространство Гёльдера), а плоская волна распространяется под углом $\pi - {{\theta }_{0}}$ по отношению к оси симметрии OZ. Тогда математическая постановка задачи рассеяния в рамках ОНО может быть записана в следующем виде, включая систему уравнений Максвелла для внешнего и внутреннего полей

граничные условия сопряжения для полей на поверхностях частиц следующего вида: и условий излучения на бесконечности для рассеянного поля Здесь $\{ {{{\mathbf{E}}}_{e}},\;{{{\mathbf{H}}}_{e}}\} $ – рассеянное поле в ${{D}_{e}}$, $\{ {{{\mathbf{E}}}_{i}},\;{{{\mathbf{H}}}_{i}}\} ,\;i = 1,2$ – полные поля внутри каждой частицы, ${{{\mathbf{n}}}_{{1,2}}}$ – единичные нормали к поверхностям $\partial {{D}_{{1,2}}} \in {{C}^{{(2,\alpha )}}}$, кроме того, внутренние электрические поля состоят из поперечных и продольных полей, то есть ${{{\mathbf{E}}}_{i}} = {\mathbf{E}}_{i}^{T} + {\mathbf{E}}_{i}^{L},$ $\operatorname{div} {\mathbf{E}}_{i}^{T} = 0,$ $\operatorname{rot} {\mathbf{E}}_{i}^{L} = 0,$ $k = {\omega \mathord{\left/ {\vphantom {\omega c}} \right. \kern-0em} c}$, а характеристики среды выбраны таким образом, что $\operatorname{Im} {{\varepsilon }_{e}},{{\mu }_{e}} = 0$, $\operatorname{Im} {{\varepsilon }_{i}},{{\mu }_{i}} \leqslant 0$, $\operatorname{Im} {{\varepsilon }_{L}} \leqslant 0$. При этом предполагается, что временная зависимость выбрана в виде $\exp \left\{ {j\omega t} \right\}$.

Параметры $\eta $ и ${{\varepsilon }_{L}}$ описывают продольную составляющую внутреннего поля ${\mathbf{E}}_{i}^{L}$ и будут определены в соответствии с моделью ОНО ниже. Следует отметить, что из формулировки задачи (1) непосредственно вытекает, что продольная компонента поля, во-первых, локализована строго внутри частицы, и, во-вторых, не вносит вклад в магнитное поле ${{{\mathbf{H}}}_{i}}$, так как $\operatorname{rot} (\operatorname{grad} \Psi ) = 0$. Введение в рассмотрение дополнительно продольного поля ${\mathbf{E}}_{i}^{L}$, подлежащего определению, требует наличия дополнительного граничного условия. Оно формулируется как обращение в нуль нормальной компоненты тока проводимости $\left( {{\mathbf{J}} \cdot {\mathbf{n}}} \right) = 0$ на поверхностях $\partial {{D}_{{1,2}}}$. В силу непрерывности тангенциальных компонент полного поля ${\mathbf{H}}$ (1) на поверхностях $\partial {{D}_{{1,2}}}$ легко показать, что в этом случае нормальная компонента ${\mathbf{n}} \cdot \operatorname{rot} {\mathbf{H}} = 0$ также непрерывна, отсюда из первого уравнения Максвелла (1) получается искомый скачок нормальной компоненты поля ${\mathbf{E}}$, который и выступает в качестве дополнительного граничного условия (1) на $\partial {{D}_{{1,2}}}$ [23].

3. ВЫЧИСЛИТЕЛЬНАЯ СХЕМА МЕТОДА ДИСКРЕТНЫХ ИСТОЧНИКОВ

Для решения задачи дифракции (1) будем использовать МДИ [25, 26]. Он обладает определенными преимуществами перед другими методами: МДИ не требует ни генерации сеток, ни использования процедур интегрирования по поверхности рассеивателя. Позволяет в аналитическом виде получать как ближнее, так и дальнее поля, причем решая задачу одновременно для всего набора углов падения плоской волны и поляризаций. Позволяет проводить анализ рассеивателей, обладающих большими значениями диэлектрической проницаемости, посредством выбора различного числа дискретных источников для представления внешнего и внутреннего полей, что оказывается весьма существенным применительно к данному случаю. Отличительной особенностью метода является то обстоятельство, что он дает возможность получать апостериорную оценку погрешности полученного решения, что позволяет контролировать реальную сходимость приближенного решения к точному.

В нашем случае поле P/S поляризованной волны может быть записано в виде

здесь ${{n}_{e}} = \sqrt {{{\varepsilon }_{е }}{{\mu }_{e}}} ,\;{{k}_{e}} = k{{n}_{e}}$.

При построении приближенного решения на основе МДИ будем существенно учитывать осевую симметрию и поляризацию плоской волны [25]. В основу представления для внешнего $\{ {{{\mathbf{E}}}_{e}},\;{{{\mathbf{H}}}_{e}}\} $ и внутреннего поперечного поля $\left\{ {{\mathbf{E}}_{i}^{T},\;{{{\mathbf{H}}}_{i}}} \right\}:\operatorname{div} {\mathbf{E}}_{i}^{T} = 0$ положим следующие векторные потенциалы

здесь $Y_{m}^{e}\left( {\xi ,z_{n}^{e}} \right) = h_{m}^{{\left( {2{\text{ }}} \right)}}\left( {{\text{ }}{{k}_{e}}{{R}_{{\xi {\text{ z}}_{n}^{e}}}}} \right){{\left( {{\rho \mathord{\left/ {\vphantom {\rho {{{R}_{{\xi {\text{ z}}_{n}^{e}}}}}}} \right. \kern-0em} {{{R}_{{\xi {\text{ z}}_{n}^{e}}}}}}} \right)}^{m}},$ $Y_{m}^{i}\left( {\xi ,z_{n}^{i}} \right) = {{j}_{m}}\left( {{\text{ }}{{k}_{i}}{{R}_{{\xi {\text{ z}}_{n}^{i}}}}} \right){{\left( {{\rho \mathord{\left/ {\vphantom {\rho {{{R}_{{\xi {\text{ z}}_{n}^{i}}}}}}} \right. \kern-0em} {{{R}_{{\xi {\text{ z}}_{n}^{i}}}}}}} \right)}^{m}}$, $\;\xi = (\rho ,z){\text{ }}$ – точка в полуплоскости $\varphi = {\text{const}}$, $R_{{\xi {\text{ }}z_{n}^{{e,i}}}}^{2} = {{\rho }^{2}} + {{\left( {z - z_{n}^{{e,i}}} \right)}^{2}}$, $\;h_{m}^{{\left( 2 \right)}}$ – сферическая функция Ханкеля, $\;{{j}_{m}}$ – сферическая функция Бесселя, ${{k}_{{e,i}}} = k\sqrt {{{\varepsilon }_{{e,i}}}{{\mu }_{{e,i}}}} $, $\left\{ {z_{n}^{{e,i}}} \right\}_{{n = 1}}^{{{{N}_{{e,i}}}}}$ – координаты дискретных источников, расположенных внутри частицы на оси 0z.

Для построения приближенного решения необходимо определить значения диэлектрической проницаемости ${{\varepsilon }_{L}}$ и продольного волнового числа ${{k}_{L}}$. В рамках ОНО ${{\varepsilon }_{L}}$ определяется как ${{\varepsilon }_{L}} = {{\varepsilon }_{i}}(\omega ) - {{\omega _{p}^{2}} \mathord{\left/ {\vphantom {{\omega _{p}^{2}} {\left( {j\gamma \omega - {{\omega }^{2}}} \right)}}} \right. \kern-0em} {\left( {j\gamma \omega - {{\omega }^{2}}} \right)}}$, где ${{\varepsilon }_{i}}(\omega )$ – диэлектрическая проницаемость металла с учетом частотной дисперсии, ${{\omega }_{p}}$ – плазменная частота для данного металла, $\gamma $ – коэффициент затухания в среде [17]. Скалярный потенциал для продольного поля удовлетворяет следующему уравнению Гельмгольца: $\left( {{{\nabla }^{2}} + k_{L}^{2}} \right)\Psi (M) = 0$, а величина продольного волнового числа определяется через значение $\eta $, как $k_{L}^{2} = {{\varepsilon }_{i}}(\omega ){\text{/}}{{\eta }^{2}}$. При этом в рамках ОНО $\eta $ представляется в следующем виде ${{\eta }^{2}} = {{{{\varepsilon }_{L}}\left( {{{\beta }^{2}} + D\left( {\gamma + j\omega } \right)} \right)} \mathord{\left/ {\vphantom {{{{\varepsilon }_{L}}\left( {{{\beta }^{2}} + D\left( {\gamma + j\omega } \right)} \right)} {\left( {{{\omega }^{2}} - j\gamma \omega } \right)}}} \right. \kern-0em} {\left( {{{\omega }^{2}} - j\gamma \omega } \right)}}$. Здесь β – гидродинамическая скорость в плазме, связанная со скоростью Ферми ${{v}_{F}}$ соотношением ${{\beta }^{2}} = 3{\text{/}}5v_{F}^{2}$, D – коэффициент диффузии электронов [23]. Отметим, что при $\beta \to 0,\;D \to 0$ поле внутри частицы становится чисто поперечным, так как ${{k}_{L}} \to \infty {\text{ }} \Rightarrow \Psi \to 0$.

Для случая Р-поляризации продольное поле строится на основе следующих скалярных потенциалов, представляющих собой частные решения уравнения Гельмгольца с нелокальным волновым числом $\left( {{{\nabla }^{2}} + k_{L}^{2}} \right){{\Psi }_{{mn}}}(M) = 0$

Итак, приближенное решение в рамках МДИ для Р-поляризованной плоской волны (2) может быть представлено в следующем виде:

В случае S-поляризации скалярный потенциал для продольного поля будет иметь вид

тогда приближенное решение для S-поляризации может быть записано, как Сравнивая представления (5) и (7), замечаем отсутствие вклада продольного поля ${\mathbf{E}}_{i}^{{LN}}$ в независящую от φ гармонику. Это является следствием того обстоятельства, что в этом случае присутствует лишь азимутальная компонента электрического поля ${{E}_{\varphi }}$ и отсутствует его нормальная компонента.

Подставляя представления для полей (5), (7) в (1), легко убедиться, что они аналитически удовлетворяют всем условиям граничной задачи (1), за исключением условий сопряжения на поверхностях частиц $\partial {{D}_{{1,2}}}$. Используя условия сопряжения для полей, получаем для определения неизвестных амплитуд ДИ следующие соотношения:

Схема вычислительного алгоритма определения амплитуд ДИ $\left\{ {p_{{mn}}^{{e,i}};p_{{mn}}^{{L,i}};q_{{mn}}^{{e,i}};r_{n}^{{e,i}};r_{n}^{{L,i}}} \right\}$ заключается в использовании обобщенного метода коллокаций [27]. Точки коллокаций $\left\{ {{{\xi }_{l}}} \right\}_{{l = 1}}^{K}$ равномерно покрывают образующие поверхностей частиц, а число источников выбирается так, чтобы линейная система для каждой Фурье гармоники оказывалась переопределенной $5K > 2\left( {{{N}_{e}} + {{N}_{i}}} \right) + {{N}_{L}}$. Далее, проводится QR-факторизация матриц и последовательное псевдорешение систем для каждой гармоники для всего набора углов падения волны ${{\theta }_{0}}$ [28]. При необходимости вычисляется невязка граничных условий в промежуточных точках по отношению к точкам коллокаций [25].

Определив амплитуды ДИ, легко вычислить θ, φ компоненты диаграммы направленности рассеянного поля ${\mathbf{F}}(\theta ,\varphi )$ на единичной сфере $\Omega = \left\{ {0 \leqslant \theta \leqslant \pi ;\;0 \leqslant \varphi \leqslant 2\pi } \right\}$

Для случая P-поляризации они принимают следующий вид: Совершенно аналогично компоненты диаграммы для случая S-поляризации будут Как следует из (9), (10), компоненты диаграммы направленности представляют собой комбинации элементарных функций, что существенно облегчает вычисления.

4. ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ

Нас будет интересовать сечение рассеяния, которое представляет собой суммарную интенсивность рассеянного поля

здесь DSC на единичной сфере определяется как Размерность DSC – μм². Кроме того, мы будем определять сечение экстинкции [29], которое применительно к нашему случаю (9), (10) принимает вид В качестве вещества частиц, расположенных в воде ($\sqrt {{{\varepsilon }_{e}}} = 1.33$), будем рассматривать золото (Au), для которого соответствующие параметры, необходимые для определения величин ${{k}_{L}}$ и ${{\varepsilon }_{L}}$, имеют значение $\hbar {{\omega }_{p}} = 9.02\;{\text{eV}},$ $\hbar \gamma = 0.071\;{\text{eV}},$ ${{v}_{F}} = 1.39\;\mu m{\text{/}}\sec ,$ $D = 1.90 \times {{10}^{8}}\;\mu {{m}^{2}}{\text{/}}\sec $ [23]. Задавая длину волны внешнего возбуждения λ, вычисляя соответствующее значение $\omega $, легко определить значения нелокальных параметров ${{k}_{L}}$ и ${{\varepsilon }_{L}}$ по приведенным выше формулам. Величина диэлектрической проницаемости для золота ${{\varepsilon }_{i}}(\omega )$ определялась с учетом частотной дисперсии металла [30]. Следует также отметить, что несмотря на имеющиеся публикации, связанные с анализом рассеивающих свойств парных частиц, в большинстве из них рассматриваются либо цилиндрические структуры (2D), либо сферические составляющие частицы, а также отсутствует исследование влияния деформации на характеристики рассеяния и ближнего поля [31]–[34].

В качестве формы частиц будем рассматривать вытянутые сфероиды, эквиобъемный диаметр каждого D = 12 нм, задавая в качестве исходных параметров соотношение осей сфероидов $r = a{\text{/}}b > 1$ и величину расстояния между частицами – d, т.е. величину зазора. Прежде чем переходить к описанию численных результатов, отметим, что основной вычислительной проблемой при учете продольных волн является весьма существенное различие величин продольного ${{k}_{L}}$ и поперечного волновых ${{k}_{e}}$ чисел. Например, в рассматриваемом ниже диапазоне длин волн это различие доходило до значения ${{\left| {{{k}_{L}}} \right|} \mathord{\left/ {\vphantom {{\left| {{{k}_{L}}} \right|} {{{k}_{e}}}}} \right. \kern-0em} {{{k}_{e}}}} \approx 200$, что определяет существенно различный характер поведения полей. В то время, как для поперечного волнового числа отношение ${{\left| {{{k}_{i}}} \right|} \mathord{\left/ {\vphantom {{\left| {{{k}_{i}}} \right|} {{{k}_{e}}}}} \right. \kern-0em} {{{k}_{e}}}}$на превышало 20. Таким образом, оказывается, что продольная волна осциллирует в 10 раз чаще, чем поперечная, что влечет за собой существенное возрастание необходимого для достижения приемлемой точности числа точек коллокаций K по сравнению с локальным случаем [25].

В нашей предыдущей публикации [35], посвященной классической модели линейного кластера плазмонных частиц, было установлено, что можно управлять положением и амплитудой ПР в оптическом диапазоне, варьируя три основных параметра: D – эквиобъемный диаметр частиц, r – вытянутость, d – величину зазора. За счет этого можно обеспечить нужное расположение ПР в широком диапазоне длин волн.

Перейдем к описанию результатов моделирования. На фиг. 1 приведены спектральные оптические характеристики: сечение экстинкции $\sigma _{{ext}}^{P}{\text{ }}$ и рассеяния $\sigma _{{scs}}^{P}$ для сфероидов с соотношением осей $r = 2$ и зазором d = 1 нм при направлении падения плоской волны перпендикулярно оси сфероидов (${{\theta }_{0}} = 90^\circ $). Из рисунка видно, что максимумы ПР для GNOR сдвинуты в коротковолновую область (blue shift) по сравнению с локальным случаем (LRA). Кроме того, очевидно существенное увеличение различия в амплитудах ПР между сечением экстинкции и сечением рассеяния для GNOR по сравнению с локальным случаем.

Фиг. 2 посвящена анализу сечения рассеяния (11) для сфероидов с $r = 2$ при различных значениях размеров зазора d. Как видно из рисунка, при уменьшении зазора различие между локальным и нелокальным ПР усиливается, при этом возрастают как величина сдвига максимума, так и различие в амплитудах. На фиг. 3 приведены значения $\sigma _{{scs}}^{P}$ при зазоре d = 1 нм для различных значений вытянутости сфероидов $r$. Из рисунка явствует, что увеличение вытянутости ведет к возрастанию амплитуды ПР. Таким образом, следует заключить, что альтернативным уменьшению зазора может служить путь по увеличению вытянутости составляющих элементов, который может привести к лучшим результатам в смысле возрастания интенсивности рассеяния. Аналогичные предыдущему результаты, но для относительной интенсивности поля ${{\left| {{{{\mathbf{E}}}_{e}} + {{{\mathbf{E}}}_{0}}} \right|}^{2}}{\text{/}}{{\left| {{{{\mathbf{E}}}_{0}}} \right|}^{2}}$ в центре зазора между частицами, можно видеть на фиг. 4, 5. Здесь также видна альтернатива уменьшению зазора для усиления ближнего поля (т.н. hot spot), которая состоит в увеличении вытянутости сфероидов.

5. ЗАКЛЮЧЕНИЕ

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