Журнал вычислительной математики и математической физики, 2021, T. 61, № 4, стр. 580-607

ВВЕДЕНИЕ

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

Оптические свойства поверхностных плазмонов очень успешно описываются на основе системы уравнений Максвелла. Это описание предполагает классическое поведение как электронов, так и электромагнитных полей в металлических наноструктурах. Наиболее привлекательная особенность поведения наноплазмонных частиц состоит в том, что частицы меньшего размера могут реализовывать усиление поля большее, чем частицы большего диаметра (см. [7]). По мере развития технологического прогресса размер элементов плазмонных структур неуклонно уменьшался. Это привело к тому, что классического описания поведения электромагнитного поля оказалось недостаточно, так как стали проявляться такие квантово-механические явления, как нелокальность и туннельные эффекты, а также двойственность волна–частица. Учет этих эффектов создает много новых возможностей для расширения границ фундаментальной науки и прикладных квантовых технологий (см. [8]).

Металлодиэлектрические структуры привлекают все более пристальное внимание, благодаря совершенствованию технологии их синтеза, возникшей в конце XX века (см. [9]). Такие структуры обеспечивают большую гибкость в управлении плазмонным резонансом (ПР), а также в достижении желаемого усиления электромагнитных полей. Это происходит в результате появления дополнительных параметров по сравнению с однородными частицами. Действительно, в случае однородных частиц манипулирование положением ПР в частотной области и его амплитудой достигается за счет размеров, материала и формы частиц. В то время, как при использовании гибридных частиц появляется дополнительный слой, который существенно расширяет возможности управления. К подобным структурам относятся, в частности, многослойные сферические и сфероидальные оболочки (см. [10]), которые в последние годы представляют повышенный интерес для исследователей в области наноплазмоники. В частности, частицы ядро–оболочка, построенные из композитных наноматериалов, стали весьма востребованными для катализа (см. [11]), сенсоров (см. [12]), фотометрического усиления (см. [13]) и солнечных элементов (см. [14]).

Фундаментальной научной проблемой в рамках квантовой плазмоники является проблема разработки и реализации наноразмерных источников когерентного излучения. Идея состоит в том, чтобы использовать плазмонные поля вместо фотонных, применяемых в обычных лазерах. Дело в том, что плазмонные поля позволяют преодолеть дифракционное ограничение размера лазера. Плазмонные нанолазеры, основанные на использовании слоистых 3D резонаторов, имеют преимущества наноразмера, низкого энергетического порога и супермалого времени отклика (см. [15]). Плазмонный нанолазер (ПН) носит название SPASER (Surface Plasmon Amplification by Simulated Emission of Radiation, далее – “спасер”) (см. [16], [17]). Первыми, кто предложил концепцию лазерного резонатора, основанного на поверхностных плазмонах, взаимодействующих с усиливающей средой, были Сударкин и Демкович (см. [18]). Теория спасера была впервые разработана Стокманом и Бергманом в 2003 г. (см. [19]). Первая демонстрация ПН была реализована Ногиновым и соавт. (см. [20]). Группа Ногинова в действующем прототипе спазера использовала в качестве резонатора одиночную золотую наночастицу сферической формы диаметром 14 нм, заключенную в кварцевую оболочку и расположенную на поверхности стеклянной призмы. Одна из возможных конструкций спасера состоит из наночастиц благородного металла, выступающих в роли нанорезонаторов, заключенных в усиливающую среду (см. [15]). В настоящее время продолжаются многочисленные исследования и разработки различных перспективных схем плазмонных нанолазеров (см. [17], [21]).

При изучении оптических свойств металлодиэлектрических плазмонных наноструктур основное внимание уделялось теоретическому анализу их спектральных свойств в дальнем поле, особенно таких характеристик, как рассеяние и экстинкция. Известно, что из-за небольшого размера этих структур по сравнению с длиной оптической волны дипольное приближение для полей представляется достаточным при описании полей в дальней зоне, что делает теоретическое исследование относительно простым. Однако существуют некоторые важные приложения, такие как молекулярная флуоресценция (см. [22]), химический катализ (см. [11]), фототермическое усиление (см. [23]) и 3D резонатор ПН (см. [17]), которые существенно определяются взаимодействием с ближним полем плазмонных структур. В этих случаях требуются более продвинутые теоретические методы для моделирования происходящих физических процессов, особенно в ближнем поле.

В рамках классического описания полей предполагается, что смещение в каждой точке зависит от приложенного электрического поля в той же точке. В этом случае среда описывается диэлектрической проницаемостью, как функцией, зависящей только от частоты возбуждающего поля. Дело в том, что когда характерный размер металлических наночастиц становится сравнимым с длиной волны Ферми электронов в этом металле (~5 нм для золота и серебра), возникает так называемая пространственная нелокальность металла (см. [24]). В этом случае электрическое поле нелокально зависит от плотности тока, индуцированного внешним полем, всюду внутри частицы. Последнее обстоятельство ведет к тому, что для частиц размером менее 10 нм требуется полное квантово-механическое описание. Например, Time-dependent density functional theory (TDDFT) или теория функционала плотности, зависящей от времени, которая описывает коллективное движение электронов, моделируя поведение каждого электрона на основе уравнения Шрёдингера. Она хорошо подходит для объяснения экспериментальных результатов, но применима лишь к плазмонным частицам диаметром всего в несколько нанометров (см. [25]). Суть в том, что когда размер плазмонных структур становится меньше длины свободного пробега возбужденных электронов, столкновениями между электронами нельзя пренебрегать. Следовательно, движение проводящих электронов будет связано не только с полем, приложенным в данной точке, но также и с полями в других точках. Таким образом, проявляется эффект нелокальности (ЭН). С целью описания оптических свойств плазмонных наночастиц были разработаны различные методы, учитывающие квантовые эффекты, встроенные в рамки классической электромагнитной теории Максвелла. Они также известны как квазиклассические подходы. Одной из наиболее популярных моделей, рассматривающих подобные нелокальные эффекты, является гидродинамическая модель Друде (ГМД) (см. [26], [27]). ГМД описывает поведение электронов внутри металла как жидкость, подчиняющуюся законам гидродинамики. Подобное рассмотрение приводит к необходимости учета возникающего продольного электрического поля, которого нет в классической модели Максвелла ($\operatorname{div} {\mathbf{E}}$ = 0).

В рамках ГМД плазмонное затухание обычно моделируется зависящей от геометрии и материала частицы диэлектрической проницаемостью (см. [28]), в которой частота затухания корректируется в соответствии с уменьшенной средней длиной свободного пробега электронов. Однако последнее обстоятельство не позволяет непосредственно использовать ГМД для исследования частиц несферической формы. Подобное неудобство было преодолено в рамках разработанной модели Обобщенного нелокального отклика (ОНО) (см. [29]). В рамках модели ОНО зависящее от размера затухание плазмона возникает естественным образом, благодаря дополнительной составляющей гидродинамического описания индуцированных зарядов в металле, а именно, коэффициенту диффузии электронов. Модель ОНО в последнее время достигла значительных успехов в области нанофотоники и наноплазмоники (см. [30]). Эта теория оказывается особенно эффективной для одновременного описания, как эффекта нелокального экранирования, так и затухания Ландау (см. [31]). В настоящее время квазиклассические модели квантовых эффектов в наноплазмонике являются наиболее востребованными, так как позволяют правильно описывать поведение оптических характеристик плазмонных наночастиц, не прибегая к чисто квантовым подходам, численные реализации которых обладают исключительно высокими требованиями к временным и компьютерным ресурсам.

Чтобы обеспечить надежные, точные и контролируемые результаты моделирования, численная электромагнитная модель должна учитывать все особенности взаимодействия в плазмонных структурах вплоть до измерения Ангстрема. Обзор различных численных методов анализа наноплазмонных структур можно найти в [32]–[34]. Среди наиболее популярных подходов следует отметить следующие:

1. Прямые методы, которые непосредственно применяются к системе уравнений Максвелла. В эту категорию включаются метод конечных разностей во временной области (FDTD) (см. [35]) и метод конечных элементов (FEM) (см. [36]), работающий в частотной области. Метод FDTD – один из самых популярных методов в нанофотонике из-за простоты его концепции и способности охватывать широкий круг задач нанофотоники. Сравнительно недавно была предложена модификация этого метода – разрывный метод Галеркина во временной области (DGTD). Он обладает несколькими привлекательными особенностями по сравнению с FDTD, такими как простота адаптации к сложной геометрии рассеивателя и составу материалов, более высокая точность и естественный параллелизм (см. [37]). FEM – еще один популярный метод, действующий в частотной области, который позволяет точно вычислять электромагнитные поля. В последние годы был предложен гибридный разрывный метод Галеркина (HDG) для решения задач в частотной области. Он вобрал в себя практически все преимущества разрывных методов Галеркина, вместе с тем, по вычислительной сложности он приближается к FEM (см. [38], [9]). Следует отметить, что перечисленные методы требуют дополнительных усилий для вычисления полей в дальней зоне.

2. Полуаналитические объемные методы: приближение дискретными диполями (DDA) (см. [34]) и метод объемных интегральных уравнений (VIE) (см. [33]). Данные методы также требуют дополнительных преобразований для вычисления полей в дальней зоне.

3. Полуаналитические поверхностные методы: метод поверхностных интегральных уравнений (SIE) (см. [33]), метод T-матриц (см. [40]), метод множественных мультиполей (MMP) (см. [41]) и метод дискретных источников (см. [42]).

В настоящее время имеется ряд статей, в которых описывается численная схема для анализа наноплазмонных структур с эффектом нелокальности (см. [43]–[46]). В то же время большинство численных результатов получены лишь для двумерных конструкций, т.е. для случая задач рассеяния, обладающего цилиндрической симметрией. Исключение составляют работы, связанные с исследованием частиц, в том числе гибридных, обладающих сферической симметрией (см. [30]). Что касается произвольных геометрий рассеивателей, то имеющиеся численные результаты в большинстве получены с использованием пакета COMSOL Multiphysics (см. [47]).

Для плазмонных структур, состоящих из однородных или слоистых плазмонных частиц, поверхностные методы представляются наиболее подходящими (см. [33]). Причина в том, что включение нелокального отклика требует решения векторного уравнения Гельмгольца с очень большим волновым числом, которое может превосходить классическое на два порядка. Большинство прямых методов ограничиваются дискретизацией шага порядка 0.2 нм, что недостаточно для моделирования нелокального отклика. Фактически дискретизация определяется длиной волны Ферми внутри плазмонного материала, что требует дискретизации шага менее 1 Å для достижения разумной точности (см. [48]). Кроме того, поверхностные методы позволяют решать задачи рассеяния сразу для всего набора внешних возбуждений и не требуют дополнительных преобразований из ближней зоны в дальнюю для вычисления рассеянных полей.

2. КОНЦЕПЦИЯ МЕТОДА ДИСКРЕТНЫХ ИСТОЧНИКОВ

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

с граничным условием на поверхности локального рассеивателя ${{D}_{i}}$:

Здесь $\hat {L}~$ – матричный оператор граничной задачи рассеяния (система уравнений Гельмгольца в случае акустических волн, система Максвелла в электромагнитном случае, система уравнений упругости для задач сейсмики и пр.), $\hat {l}$ – матричный оператор граничных условий (условий сопряжения), $\partial {{D}_{i}} \in {{C}^{{(1,\alpha )}}}$ – замкнутая гельдерова поверхность, ограничивающая односвязную область D_i и U⁰ – возбуждающее поле.

При рассмотрении граничных задач теории дифракции существенно, что рассеивающий объект расположен на значительном удалении от источников первичного поля. В этом случае возбуждающее поле на поверхности рассеивателя представляет собой аналитическую функцию вещественных переменных. Это обстоятельство позволяет использовать концепцию квазирешения для построения приближенного решения (см. [49]).

Суть концепции квазирешения в рассматриваемом случае заключается в построении аналитической конструкции для приближенного решения граничной задачи (2.1)–(2.3) ${{U}^{\delta }}$, которое должно удовлетворять аналитически уравнению (2.1) и условию на бесконечности (2.3), а граничное же условие на поверхности рассеивателя удовлетворяется приближенно, как

Последнее обстоятельство достигается с использованием соответствующего вычислительного алгоритма, позволяющего численно определять параметры выбранного представления для приближенного решения так, чтобы обеспечить заданную величину невязки на поверхности рассеивателя (2.4). Таким образом, граничная задача (2.1)–(2.3), первоначально сформулированная во всем пространстве, сводится к задаче аппроксимации полей на поверхности препятствия. Фундаментальную роль при решении задачи аппроксимации играет то обстоятельство, что выполнение условия (2.4) при достаточно малом δ обеспечивает близость в равномерной метрике приближенного решения к точному всюду вне рассеивателя.

Действительно, в предположении существования единственного решения граничной задачи дифракции (2.1)–(2.3) существует соответствующий тензор Грина $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\leftrightarrow}$}} {G} (M,{{M}_{0}})$, удовлетворяющий условию

и условиям не бесконечности. Откуда, используя формулу Грина, нетрудно получить соотношение корректности для граничной задачи рассеяния вида где d – любой компакт в D_e. Итак, чтобы обеспечить сходимость приближенного решения к точному решению в равномерной метрике вне рассеивателя, достаточно аппроксимировать с произвольной точностью граничное условие (2.2) в метрике L₂. Таким образом, решение исходной граничной задачи дифракции сводится к задаче аппроксимации граничных значений внешнего возбуждения на поверхности локального препятствия D_i.

Вопрос об обеспечении выполнения условия (2.4) для произвольного $\delta \to 0~$ является центральным при построении приближенного решения. Это может быть реализовано различными способами. Большую роль в этом плане играет носитель множества дискретных источников. В [50], [51] подробно изложена схема построения полных и замкнутых систем МДИ, которые используются при решении задач рассеяния. Применительно к системе Максвелла в качестве носителя может быть выбрана замкнутая поверхность, тогда в качестве источников достаточно выбрать три ортогональных электрических диполя, распределенных по поверхности. В случае разомкнутой поверхности достаточной оказывается система тангенциальных электрических и магнитных диполей. При выборе в качестве носителя отрезка прямой “базисная” система функций будет представлять собой систему распределенных мультиполей низшего порядка (Lowest Order Distributed Multipoles) (см. [52]).

Предположим, что выбран какой-то носитель дискретных источников γ. Это может быть поверхность или ее часть, часть кривой в пространстве или отрезок прямой. Тогда представление для приближенного решения задачи (2.1)–(2.3) можно записать как

где ${{\hat {T}}_{\gamma }}$ – некоторый линейный оператор, определенный на ${{L}_{2}}(\gamma )$. Нас будет интересовать его свойства в совокупности с граничным оператором (2.2) на поверхности $\partial {{D}_{i}}$, т.е. свойства оператора Поскольку ${{U}^{\delta }}(M)$ является аналитической функцией всюду вне носителя, то достаточно будет показать, что замыкание области значений оператора $\hat {l}{{\hat {T}}_{\gamma }}$ совпадает с пространством ${{L}_{2}}(\partial {{D}_{i}})$. Для этого достаточно показать, что ядро сопряженного к $\hat {l}{{\hat {T}}_{\gamma }}~$ оператора пусто (см. [53]). Обоснование последнего факта обычно достигается сведением к однородным уравнениям метода нулевого поля, которые, в свою очередь, эквивалентны однородной граничной задаче рассеяния (2.1)–(2.3), имеющей нулевое решение (см. [52], [54]).

Теоретические основы метода вспомогательных источников (Method of Auxiliary Sources – MAS) были заложены в работах В.Д. Купрадзе (см. [55]) и М.А. Алексидзе (см. [56]). Различные модификации этого подхода использовались в многочисленных практических приложениях, как в России, так и за рубежом. В приложении к двумерным задачам следует отметить работы Р. Заридзе [57], [58], K. Яшуура (K. Yasuura) [59], [60], И. Левиатана (Y. Leviatan) [61], [62] и Г. Фикиориса (G. Fikioris) [63]. В случае пространственных задач рассеяния, не обладающих осевой симметрией, этот подход более известен как метод множественных мультиполей (Multiple Multipole Method – MMP). Наряду с работами авторов [64], [65] здесь необходимо отметить школу Х. Хафнера (Ch. Hafner) (см. [66], [67]), который и ввел понятие MMP. С последними инновациями в области MAS–MMP можно ознакомиться в [68].

МДИ представляется одним из наиболее эффективных и гибких инструментов для построения квазирешения в задачах рассеяния волн (среди прочих численно-аналитических методов). В рамках МДИ приближенное решение строится как конечная линейная комбинация полей диполей и мультиполей с неизвестными амплитудами. Это представление аналитически удовлетворяет уравнению (2.1) и условию излучения на бесконечности. Неизвестные амплитуды дискретных источников (ДИ) должны определяться из граничного условия (2.2). Одна из наиболее привлекательных особенностей МДИ заключается в гибком выборе полей ДИ, которые используются для построения приближенного решения. В рамках МДИ нет никаких ограничений на выбор источников поля. Поля ДИ должны лишь удовлетворять уравнению (2.1), условию излучения и образовывать полную и замкнутую функциональную систему на поверхности локального препятствия. В этом случае обеспечивается сходимость приближенного решения к точному решению, т.е. выполнение (2.4) для произвольного значения δ. В силу (2.5) возможно контролировать погрешность приближенного решения, оценивая невязку поля на поверхности рассеивателя. Это обстоятельство приобретает особое значение в случае необходимости вычисления величины полного поля вблизи поверхности рассеивателя.

В случае рассмотрения осесимметричных структур наибольшее распространение получила система распределенных мультиполей низшего порядка, впервые реализованная в [69]. В этом случае источники как внешнего, так и внутреннего поля располагаются на отрезке оси симметрии внутри рассеивателя. Подобная схема позволяет эффективно исследовать вытянутые тела вращения (см. [70]). Причем для представления внутреннего поля используются функции Бесселя. В случае сплюснутых рассеивателей была предложена процедура продолжения полей с оси в комплексную плоскость по координате источника (см. [71]). Расположение источников в комплексной плоскости позволило рассматривать сплюснутые рассеиватели с большим соотношением осей (см. [72], [73]). Кроме того, при рассмотрении осесимметричных структур решение строится с учетом поляризации исходного излучения. Последнее означает, что оно имеет различный вид для P и S поляризаций. Вместе с тем численная схема определения амплитуд ДИ организована таким образом, чтобы решать единую систему уравнений для P и S поляризаций, кроме гармоник, не зависящих от азимутальной переменной (см. [74]).

Большое число работ по МДИ связано с рассмотрением осесимметричных структур в присутствии слоистой среды (см. [75]–[77]). В этом случае представление для решения в области вне локального рассеивателя строится на основе тензора Грина слоистой среды (см. [78]). Данное обстоятельство позволяет сохранить численную схему определения амплитуд ДИ без изменений. Интересно, что подобный подход был предложен еще в [79], но в те годы оказался невостребованным.

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

3. ТЕОРИЯ НЕЛОКАЛЬНОГО ОТКЛИКА

В 70-х годах прошлого столетия было замечено, что наноразмерные металлические частицы при облучении светом проявляют аномальные рассеивающие свойства. Возникла гипотеза о том, что индекс рефракции подобных частиц отличается от индекса рефракции металлических пленок, которые служат основой для измерений индексов. Пионером в изучении этого явления явился израильский ученый Р. Руппин (R. Ruppin) (см. [83]). Именно его работы заложили теоретические основы анализа влияния эффекта нелокального взаимодействия на рассеивающие свойства наноразмерных металлических частиц (см. [84], [85]). Однако в те годы потребности практического использования этих результатов не представлялись актуальными.

В данном разделе приводится краткий обзор теории нелокального отклика, следуя публикациям [24], [31]. Как известно, основным соотношением, связывающим вектор электрической индукции D и вектор напряженности электрического поля E, является соотношение

где ω – частота колебаний поля, ε₀ – диэлектрическая проницаемость вакуума и ${{\varepsilon }}({\mathbf{r}},{\mathbf{r}}{\kern 1pt} ',{{\omega }})$ – скалярная нелокальная диэлектрическая проницаемость металла. В приближении стандартного локального отклика имеет место соотношение ${{\varepsilon }}({\mathbf{r}},{\mathbf{r}}{\kern 1pt} ',{{\omega }}) = {{\delta }}({\mathbf{r}} - {\mathbf{r}}{\kern 1pt} '){{\varepsilon }}({{\omega }})$, откуда получаем где ${{\varepsilon }}({{\omega }})$ – пространственно независимая диэлектрическая проницаемость с учетом частотной дисперсии. В рамках модели нелокального отклика векторное волновое уравнение для электрического поля записывается как где $с = \frac{1}{{\sqrt {{{{{\varepsilon }}}_{0}}{{{{\mu }}}_{0}}} }}$ – скорость света в вакууме. В нелокальном случае диэлектрическая проницаемость может быть представлена в виде где $f({\mathbf{r}} - {\mathbf{r}}{\kern 1pt} ',{{\omega }})$ – скалярная функция нелокального отклика. Функция нелокального отклика предполагается симметричной и короткодействующей, откуда следует (см. [24]), что где ${{{\mathbf{r}}}^{2}} = {{x}^{2}}{{{\mathbf{e}}}_{x}} + {{y}^{2}}{{{\mathbf{e}}}_{y}} + {{z}^{2}}{{{\mathbf{e}}}_{z}}$, (${{{\mathbf{e}}}_{x}}$, $~{{{\mathbf{e}}}_{y}}$, $~{{{\mathbf{e}}}_{z}}$) – базис в декартовой системе координат, а ξ – масштаб взаимодействия нелокального отклика. Согласно этому предположению и учитывая разложение ${\mathbf{E}}({\mathbf{r}}{\kern 1pt} ',{{\omega }})$ в ряд Тейлора второго порядка в окрестности r, интеграл (3.3) может быть вычислен аналитически и векторное волновое уравнение (3.3) принимает вид Из последнего соотношения видно, что скалярный нелокальный отклик определяется Лапласианом в векторном волновом уравнении. Гидродинамическая модель Друде и модель обобщенного нелокального оптического отклика являются теориями, описывающими нелокальный отклик и приводящими к векторному волновому уравнению в форме (3.5). Это позволяет получить простые формулы для вычисления масштаба длины ξ.

В гидродинамической модели для газа свободных электронов энергия электронной плазмы выражается через электронную плотность $n({\mathbf{r}},t)$ и поле скоростей ${\mathbf{v}}({\mathbf{r}},t)$. Динамика $n({\mathbf{r}},t)$ и ${\mathbf{v}}({\mathbf{r}},t)$, обусловленная электрическим полем ${\mathbf{E}}({\mathbf{r}},t)$, получается дифференцированием энергии по этим переменным. Дифференцирование энергии по скорости и электронной плотности дает, соответственно, гидродинамическое уравнение движения и уравнение непрерывности. Эти уравнения решаются с помощью линеаризации. В результате в частотной области мы приходим к векторному волновому уравнению следующего вида (см. [31]):

и соотношению для плотности тока в виде Здесь ${\mathbf{J}}({\mathbf{r}},{{\omega }}) = - e{{n}_{0}}{\mathbf{v}}({\mathbf{r}},{{\omega }})$ – плотность тока свободных зарядов, ${{n}_{0}}$ – равновесная электронная плотность свободных электронов, ${{\beta }^{2}} = \left( {\frac{3}{5}} \right){v}_{F}^{2}$, ${{{v}}_{F}}$ – скорость Ферми, ${{\gamma }}$ – скорость затухания Друде, $\omega _{p}^{2} = \frac{{{{n}_{0}}{{e}^{2}}}}{{{{\varepsilon }_{0}}m}}$ –  плазменная частота металла, $\sigma ({{\omega }}) = - j{{\varepsilon }_{0}}\frac{{\omega _{p}^{2}}}{{\omega - j\gamma }}$ – проводимость Друде, связанная с диэлектрической проницаемостью ${{\varepsilon }}({{\omega }})$ посредством соотношения вида $\varepsilon (\omega ) = {{\varepsilon }_{b}} - j\frac{{\sigma (\omega )}}{{{{\varepsilon }_{0}}\omega }} = {{\varepsilon }_{b}} - \frac{{\omega _{p}^{2}}}{{{{\omega }^{2}} - j\gamma \omega }}$, и ${{\varepsilon }_{b}}$ – диэлектрическая проницаемость, соответствующая связанным зарядам (bound charges). Очевидно, в пределе $\beta \to 0$, т.е. ${{{v}}_{F}} \to 0,$ основное соотношение (3.7) упрощается, сводясь ко всем известному закону Ома ${\mathbf{J}} = {{\sigma }}{\mathbf{E}}$.

Гидродинамическая модель включает в себя конвективный ток, но пренебрегает диффузионными токами. Это создает проблемы при рассмотрении рассеивателей, форма которых отличается от сферической. В данном случае приходится корректировать параметр γ, описывающий скорость затухания (damping rate) электронов в металле (см. [28]). В модели обобщенного нелокального оптического отклика гидродинамическая теория дополнительно учитывает диффузию электронов. Учет электронной диффузии меняет уравнение непрерывности, которое в своей линеаризованной форме представляет собой конвекционно-диффузионное уравнение. В этой модели плотность токов свободных зарядов учитывает диффузионный член ${\mathbf{J}}({\mathbf{r}},{{\omega }}) = - e{{n}_{0}}{\mathbf{v}}({\mathbf{r}},{{\omega }}) + eD\nabla n({\mathbf{r}},{{\omega }})$, где D – коэффициент диффузии, и основное соотношение для плотности токов свободных зарядов становится (см. [29])

Для того чтобы получить векторное волновое уравнение для поля, исключим плотность тока из уравнений (3.6) и (3.8). Применяя оператор $\nabla (\operatorname{div} )$ к уравнению (3.6) и учитывая, что $(\operatorname{div} \operatorname{rot} {\mathbf{F}}) = 0$, подставим результат в (3.8), в итоге получаем

тогда уравнение (3.6) преобразуется к виду где ${{\xi }}$ – масштаб длины в обобщенной модели нелокального оптического отклика имеет вид Очевидно, что векторное волновое уравнение (3.10) такое же, как в (3.5), но в нем оператор Лапласа заменен на оператор градиента от дивергенции.

Фактически присутствие оператора Лапласа в (3.5) является следствием использования скалярной функции, описывающей нелокальный отклик $f({\mathbf{r}} - {\mathbf{r}}{\kern 1pt} ',{{\omega }})$. Из (3.3) и (3.10) следует, что в теории нелокального отклика вектор электрической индукции, включающий вклады как связанных, так и свободных зарядов имеет вид

Действительно, векторное волновое уравнение (3.10) эквивалентно уравнениям Максвелла для немагнитных сред вида Применяя преобразования $\sqrt {{{{{\varepsilon }}}_{0}}} {\mathbf{E}}({\mathbf{r}},{{\omega }}) \to {\mathbf{E}}({\mathbf{r}},{{\omega }})$ и $\sqrt {{{{{\mu }}}_{0}}} {\mathbf{H}}({\mathbf{r}},{{\omega }}) \to {\mathbf{H}}({\mathbf{r}},{{\omega }})$, можно переписать уравнения Максвелла в безразмерном виде где ${{k}_{0}} = {{\omega }}\sqrt {{{{{\varepsilon }}}_{0}}{{{{\mu }}}_{0}}} = \frac{\omega }{c}$ – волновое число в вакууме. В данном случае электрическое поле представляет собой суперпозицию вихревого поперечного и невихревого продольного полей, т.е. где $\operatorname{div} {{{\mathbf{E}}}_{{\text{T}}}}({\mathbf{r}},{{\omega }}) = 0$ и $\operatorname{rot} {{{\mathbf{E}}}_{{\text{L}}}}({\mathbf{r}},{{\omega }}) = 0$. Подставляя (3.17) в (3.10) и используя тождество $\operatorname{rot} \operatorname{rot} {\mathbf{F}} = \nabla (\operatorname{div} {\mathbf{F}}) - \Delta {\mathbf{F}}$, мы получаем, что поперечное и продольное поля удовлетворяют следующим волновым уравнениям: где поперечное и продольное волновые числа задаются соотношениями Таким образом, нелокальный отклик возбуждает дополнительно продольную волну в металлической среде. Эта волна обусловлена плотностью зарядов и не влияет на магнитное поле, так как

Присутствие продольной волны требует введения дополнительного граничного условия, наряду с обычными условиями непрерывности для тангенциальных компонент E и H на границах раздела сред. Для того чтобы получить дополнительное граничное условие, рассмотрим диэлектрическую среду 1 и металлическую среду 2 с границей раздела S. Положим, что единичный вектор нормали ${\mathbf{\hat {n}}}$ в точке r на S направлен в диэлектрическую среду 1. Вектор электрической индукции, заданный формулой (3.2), может быть записан как

где в соответствии с (3.9) нелокальный отклик порождается плотностью токов свободных зарядов J, тогда как есть вектор индукции электрического поля связанных зарядов, и P – вектор поляризации.

Из закона Гаусса для электрического поля $\operatorname{div} {{{\mathbf{D}}}_{{\mathbf{b}}}} = {{\rho }}$, где ${{\rho }}$ – плотность свободных зарядов, применяя теорему Гаусса для замкнутой поверхности, включающей границу раздела, мы получаем граничное условие ${\mathbf{\hat {n}}} \cdot ({{{\mathbf{D}}}_{{{\text{b}}1}}} - {{{\mathbf{D}}}_{{{\text{b}}2}}}) = {{{{\rho }}}_{{\text{s}}}}$, где

есть плотность свободных поверхностных зарядов. Предполагая, что плотность свободных поверхностных зарядов исчезающее мала, т.е. ${{{{\rho }}}_{{\text{s}}}} = 0$, мы получаем условие непрерывности нормальной составляющей вектора электрической индукции связанных зарядов которое даeт дополнительное граничное условие где ${{{{\varepsilon }}}_{1}}$ – диэлектрическая проницаемость диэлектрической среды 1. Если ${{{{\varepsilon }}}_{1}} \ne {{{{\varepsilon }}}_{{\text{b}}}}$, то нормальная составляющая электрического поля разрывна на границе раздела. Этот скачок электрического поля обусловлен электрическим зарядом, возникающим при поляризации связанных электронов в диэлектрике и металле. Определяя плотность связанных поверхностных зарядов по формуле ${{{{\rho }}}_{{{\text{sb}}}}} = {\mathbf{\hat {n}}} \cdot ({{{\mathbf{P}}}_{1}} - {{{\mathbf{P}}}_{2}})$, из (3.23) и непрерывности ${\mathbf{\hat {n}}} \cdot {{{\mathbf{D}}}_{{\text{b}}}}$ получаем Именно эта плотность поверхностных зарядов ответственна за многие эффекты, возникающие при учете наличия продольных полей в плазмонных наночастицах (см. [86]).

Однако из уравнения непрерывности $\operatorname{div} {\mathbf{J}} = - j{{\omega \rho }}$ мы имеем граничные условия ${\mathbf{\hat {n}}} \cdot ({{{\mathbf{J}}}_{1}} - {{{\mathbf{J}}}_{2}}) = - j{{\omega }}{{{{\rho }}}_{{\text{s}}}}$. Предполагая снова, что плотность свободных зарядов на границе раздела равна нулю, получаем условие непрерывности нормальной составляющей плотности токов свободных зарядов ${\mathbf{\hat {n}}} \cdot {{{\mathbf{J}}}_{1}} = {\mathbf{\hat {n}}} \cdot {{{\mathbf{J}}}_{2}}$. Поскольку диэлектрик не проводит ток (в диэлектриках нет свободных электронов), то ${{{\mathbf{J}}}_{1}} = 0$. Следовательно, на границе раздела диэлектрика и проводника нормальная компонента плотности тока свободных зарядов в металле обращается в нуль, т.е. ${\mathbf{\hat {n}}} \cdot {{{\mathbf{J}}}_{2}} = 0$. Другими словами, свободные электроны не выходят за поверхность металла. Условие (3.25) известно, как “hard wall boundary condition” (см. [87]). Более общие дополнительные граничные условия на поверхности между двумя плазмонными металлами можно найти в [88].

4. МЕТОД ДИСКРЕТНЫХ ИСТОЧНИКОВ: СЛОИСТАЯ НАНОЧАСТИЦА В СВОБОДНОМ ПРОСТРАНСТВЕ

Гибридные частицы, состоящие из слоев диэлектрика или полупроводника и плазмонных металлов, вызывают все больший интерес. Это связано с появлением дополнительных параметров, которые можно использовать для оптимизации плазмонного отклика. Действительно, в случае однородной частицы варьировать можно лишь ее форму, в то время как в случае слоистой структуры можно выбирать дополнительно параметры диэлектрика, его толщину и уже потом форму гибридной частицы. В классе несферических частиц “ядро–оболочка” удлиненные частицы сфероидальной формы обеспечивают простой контроль центральной частоты плазмонного резонанса (ПР) за счет изменения длины частицы (см. [89], [90]). Существует множество практических приложений, требующих оптимизации оптических свойств частиц “ядро–оболочка” с тонкой металлической оболочкой (см. [91], [92]), что обусловливает необходимость создания моделей несферических частиц, учитывающих нелокальный эффект в плазмонной металлической оболочке.

Будем рассматривать математическую модель задачи рассеяния плоской волны на осесимметричной слоистой частице с внутренней областью ${{D}_{i}}$, областью слоя ${{D}_{s}}$ и внешней неограниченной областью ${{D}_{e}}$. Пусть области разделены гладкими, $\partial {{D}_{{i,s}}} \in {{C}^{{(1,\alpha )}}}$, поверхностями с общей осью симметрии 0z. Будем полагать, что плоская волна $\{ {{{\mathbf{E}}}^{0}},{{{\mathbf{H}}}^{0}}\} $ распространяется под углом $\pi - {{\theta }_{0}}$ по отношению к оси 0z. Тогда математическая постановка задачи рассеяния в рамках ОНО может быть записана в следующем виде:

с условиями излучения на бесконечности Здесь $\{ {{{\mathbf{E}}}_{{i,s}}},{{{\mathbf{H}}}_{{i,s}}}\} $ – полные поля в областях ${{D}_{{i,s}}}$, $\{ {{{\mathbf{E}}}_{e}},{{{\mathbf{H}}}_{e}}\} $ – рассеянное поле в ${{D}_{e}}$, ${{{\mathbf{n}}}_{{i,s}}}$ – единичные внешние нормали к поверхностям $\partial {{D}_{{i,s}}}$, ${{\varepsilon }_{{i,s,e}}}$ – диэлектрические проницаемости сред в соответствующих областях, при этом $\operatorname{Im} {{\varepsilon }_{{i,e}}}$ = 0, $\operatorname{Im} {{\varepsilon }_{{s,L}}} \leqslant 0$, $k = \frac{\omega }{c}$, временная зависимость величин выбрана в виде ${\text{exp}}(j\omega t)$.

В данном случае параметры, связанные с моделью ОНО, имеют вид

Как уже отмечалось ранее, поле в оболочке ${{D}_{s}}$ разделяется на поперечную и продольную компоненты ${{{\mathbf{E}}}_{s}} = {{{\mathbf{E}}}_{{sT}}} + {{{\mathbf{E}}}_{{sL}}}$, для которых $\operatorname{div} {{{\mathbf{E}}}_{{sT}}} = 0$ и $\operatorname{rot} {{{\mathbf{E}}}_{{sL}}} = 0$. Заметим, что магнитное поле внутри металлической оболочки остается чисто поперечным.

Будем строить приближенное решение граничной задачи (1) на основе метода Дискретных источников. Разделим поле падающей плоской волны на P и S поляризации, тогда

где $\chi (x,z) = {\text{exp}}\left\{ { - j{{k}_{e}}(x\sin {{\theta }_{0}} - z\cos {{\theta }_{0}})} \right\}$, ${{k}_{e}} = k\sqrt {{{\varepsilon }_{e}}} $, а $({{e}_{x}},{{e}_{y}},{{e}_{z}})$ – декартов базис.

Будем строить поперечные поля в областях ${{D}_{{i,s,e}}}$ на основе векторных потенциалов, индуцированных источниками, распределенными вдоль оси симметрии рассеивателя, следующего вида (см. [51]):

Здесь где $h_{m}^{{(2,1)}}$ – сферические функции Ханкеля, соответствующие “уходящим” (+) и “приходящим” (–) волнам, ${{j}_{m}}$ – сферические функции Бесселя, ${{k}_{{i,s}}} = k\sqrt {{{\varepsilon }_{{i,s}}}} $, $\zeta = (\rho ,z)$, ${{\rho }^{2}} = {{x}^{2}} + {{y}^{2}}$, $R_{{{{z}_{n}}}}^{2} = {{\rho }^{2}} + {{(z - {{z}_{n}})}^{2}}$, $\{ z_{n}^{\alpha }\} _{{n = 1}}^{{N_{\alpha }^{m}}}$ – координаты дискретных источников, распределенных вдоль оси вращения 0z, число которых может различаться в зависимости от номера гармоники по φ.

В свою очередь, приближенное решение для продольных полей в оболочке строится на основе скалярных потенциалов, которые удовлетворяют уравнению Гельмгольца вида (3.19) с волновым числом $k_{L}^{2} = \frac{{{{\varepsilon }_{s}}}}{{{{\xi }^{2}}}}$ и будет иметь вид (см. [93])

Следует отметить, что внутри оболочки поля представляются в виде суммы уходящих и приходящих волн, т.е. ${{{\mathbf{E}}}_{s}} = {{{\mathbf{E}}}_{{s + }}} + {{{\mathbf{E}}}_{{s - }}}$, ${{{\mathbf{H}}}_{s}} = {{{\mathbf{H}}}_{{s + }}} + {{{\mathbf{H}}}_{{s - }}}$. Тогда конкретные представления для поперечных и продольных полей в случае поляризации P приобретают вид В случае S поляризации для построения продольных полей используются следующие потенциалы: В свою очередь, представление для приближенного решения в этом случае может быть записано как Видно, что представления для решения в различных областях ${{D}_{{i,s,e}}}~$выглядят различным образом. Особое значение имеет представление для полей внутри оболочки, которое включает в себя продольные поля следует отметить, что представления для полей (4.11), (4.13) удовлетворяют модифицированной системе уравнений Максвелла (4.1), (4.2) и условиям излучения, остается удовлетворить условиям сопряжения для полей (4.3), (4.4).

Для определения вектора неизвестных амплитуд ДИ будем использовать обобщенный метод коллокаций (см. [94]). Для этого на образующей поверхности вращения $\partial {{D}_{{i,s}}}$ выберем точки коллокаций ${{\tau }_{{\beta l}}} = ({{\rho }_{{\beta l}}},{{z}_{{\beta l}}})$, $\beta = i,s;~\;\{ {{\tau }_{{il}}}\} _{{l = 1}}^{{{{K}_{i}}}},\;\{ {{\tau }_{{sl}}}\} _{{l = 1}}^{{{{K}_{s}}}}$. И запишем условия сопряжения в виде проекционных соотношений вида

Таким образом, для каждой отдельной гармоники по φ для определения вектора неизвестных амплитуд ДИ, ${\mathbf{p}}_{m}^{N} = \{ p_{{mn}}^{{e,i,s \pm }},q_{{mn}}^{{e,i,s \pm }},r_{n}^{{e,i,s \pm }};p_{{mn}}^{{L \pm }},q_{{mn}}^{{L \pm }},r_{n}^{{L \pm }}\} ,$ размерностью $2N_{i}^{m} + 4N_{s}^{m} + 2N_{L}^{m} + 2N_{e}^{m},$ у нас имеется $5{{K}_{i}} + 5{{K}_{s}}$ уравнений. Тогда для каждой гармоники $m = \overline {0,M} $ мы имеем матрицу размерности $(5{{K}_{i}} + 5{{K}_{s}}) \times (2N_{i}^{m} + 4N_{s}^{m} + 2N_{L}^{m} + 2N_{e}^{m})$. При использовании обобщенного метода коллокаций матрица получается переопределенной (см. [51]). Сначала осуществляется регуляризация по А.Н. Тихонову в норме ${{l}_{2}}$ (см. [95]). Псевдорешение полученной системы осуществляется факторизацией QR матрицы (см. [96]) с последующим вычислением псевдорешений для всего набора внешних возбуждений. Здесь уместно заметить, что, несмотря на различие представлений для P и S поляризаций (4.11), (4.13), удается свести вычисление амплитуд ДИ к единой матрице для всех гармоник, кроме не зависящих от φ (см. [51]).

Важной характеристикой, описывающей реакцию рассеивателя на внешнее возбуждение, служит диаграмма направленности рассеянного поля. Она определяется как (см. [97])

Следуя асимптотике полей (4.11), компоненты θ, φ диаграммы направленности для P поляризации принимают следующий вид: а для поляризации S компоненты диаграммы могут быть записаны как Таким образом, определив амплитуды ДИ, можно легко вычислять компоненты диаграммы рассеяния (4.18), (4.19).

В задачах дифракции интерес представляет вычисление интенсивности рассеянного поля DSC, которая на единичной сфере определяется соотношением

а также сечение рассеяния: суммарная интенсивность рассеянного поля где Ω – единичная сфера: $~{{\Omega }} = \{ 0^\circ \leqslant \theta \leqslant 180^\circ ;\;0^\circ \leqslant \varphi \leqslant 360^\circ \} $. Размерность интенсивности и СР в силу определения (4.17) – мкм².

Если бы в основу представлений (4.18), (4.19) были положены сферические гармоники, то в силу их ортогональности на Ω СР представлялось бы, как сумма квадратов амплитуд этих гармоник. В нашем случае легко убедиться, что функции ${{\sin }^{m}}\theta exp\{ j{{k}_{e}}z_{n}^{e}\cos \theta \} $ в рамках одной гармоники Фурье для различных положений ДИ не являются ортогональными. Покажем, что тем не менее СР может быть вычислено аналитически (см. [72]). Представим СР следующим образом:

Подставляя выражения (4.18), (4.19) и интегрируя по φ, получаем, что для определения СР (4.20) достаточно уметь вычислять следующие интегралы: где ${{\gamma }_{{nl}}} = j{{k}_{e}}(z_{n}^{e} - z_{l}^{e})$. Интегрируя по частям $I_{m}^{{(0)}}$, легко убедиться, что имеет место следующее базовое рекуррентное соотношение: где ${{a}_{m}} = 4m(m - 1){\text{/}}\gamma _{{nl}}^{2}$, ${{b}_{m}} = 2m(2m - 1){\text{/}}\gamma _{{nl}}^{2}$, $I_{0}^{{\left( 0 \right)}} = [\exp ({{\gamma }_{{nl}}}) - \exp ( - {{\gamma }_{{nl}}})]{\text{/}}{{\gamma }_{{nl}}}$, $I_{1}^{{\left( 0 \right)}} = 2[\exp ({{\gamma }_{{nl}}}) + $ $ + \;\exp ( - {{\gamma }_{{nl}}})]{\text{/}}\gamma _{{nl}}^{2} - 2[\exp ({{\gamma }_{{nl}}}) - \exp ( - {{\gamma }_{{nl}}})]{\text{/}}\gamma _{{nl}}^{3}$ при ${{\gamma }_{{nl}}} \ne 0$.

Также очевидно, что $I_{m}^{{(1)}} = \frac{\partial }{{\partial \gamma }}I_{m}^{{(0)}}$, $I_{m}^{{(2)}} = \frac{\partial }{{\partial \gamma }}I_{m}^{{(1)}}$, $m \geqslant 0$. Случай ${{\gamma }_{{nl}}} = 0$ является тривиальным, так как все интегралы $I_{m}^{{(0,1,2)}}~$, $~m \geqslant 0$, берутся аналитически.

Перейдем к результатам моделирования. Будем рассматривать слоистую частицу, состоящую из диэлектрического ядра SiO₂ (n_i = 1.46) сфероидальной формы, покрытую золотой (Au) оболочкой толщиной d, располагающейся в воде (n_e = 1.333). Нас будут интересовать численные результаты анализа оптических характеристик рассеяния и интенсивности ближнего поля в области длин волн. В частности, мы будем рассматривать поведение сечения рассеяния (4.20) и сечения экстинкции

а также коэффициента усиления поля на внешней поверхности слоистой частицы $\partial {{D}_{s}}$: Напомним, что МДИ позволяет вычислять ближние поля, компоненты диаграммы рассеяния и сечения экстинкции с учетом ЭН в аналитическом виде. Соответствующие квантовые параметры для золота в рамках модели ОНО будут равняться (см. [31]) Задавая длину волны внешнего возбуждения λ и вычисляя соответствующее значение $\omega $, легко определить значения нелокальных параметров ${{\varepsilon }_{L}}$ и ${{k}_{L}}$ (4.6). Отметим, что локальное значение диэлектрической проницаемости для золота ${{\varepsilon }_{s}}$ определялось с учетом частотной дисперсии материала Au (см. [98]). Пусть эквиобъемный диаметр ядра ${{D}_{c}} = 14$ нм, а толщины оболочки выбраны равными $d = 2.5,~\;3.5$ нм.

На фиг. 1а приведены значения полного сечения рассеяния, соответствующего решению задачи дифракции плоской волны, падающей под углом ${{\theta }_{0}} = 90^\circ $ на сфероидальную частицу с соотношением осей $r = 2$, толщина оболочки при этом $d = 3.5$ нм. Волна распространяется перпендикулярно большей оси сфероида. Из фиг. 1а следует, что S и P поляризации реализуют плазмонный резонанс (ПР) в разных частях спектра, кроме того, P поляризованное излучение доминирует. Учет ЭН приводит к снижению амплитуды ПР (damping) при небольшом сдвиге в область коротких волн (blue shift). Тот же характер поведения кривых можно наблюдать на фиг. 1б, где приведены результаты для сечения экстинкции. В этом случае снижение амплитуды ПР достигает 20%.

На фиг. 2а, б изображены результаты расчетов поведения сечения рассеяния и экстинкции для сфероида в зависимости от соотношения осей при толщине пленки $d = 2.5$ нм. Как следует из рисунков, при увеличении вытянутости сфероида, т.е. увеличении объема и поверхности плазмонного слоя влияние ЭН существенно усиливается. В этом случае снижение амплитуды ПР может доходить до 2.5 раз.

На фиг. 3 представлены результаты расчета коэффициента усиления поля на внешней поверхности слоистой частицы. В этом случае учет ЭН оказывается наиболее заметным и снижение ПР может достигать значения 70%. Наконец, фиг. 4 демонстрирует зависимость коэффициента усиления от угла падения волны для значений λ, соответствующих положению ПР в области длин волн при вытянутости $r = 2$. Видно, что коэффициент монотонно возрастает и достигает своего максимума при ${{\theta }_{0}} = 90^\circ $. Аналогичные результаты имеют место для сечения экстинкции.

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

5. МЕТОД ДИСКРЕТНЫХ ИСТОЧНИКОВ: СЛОИСТАЯ НАНОЧАСТИЦА НА ПОВЕРХНОСТИ ПРОЗРАЧНОЙ ПОДЛОЖКИ

В целом ряде приложений, упомянутых выше, возникает необходимость в исследовании оптических свойств наночастиц, расположенных на диэлектрической подложке. В этом случае весьма важно строго учитывать взаимодействие между частицей и поверхностью подложки. Эти обстоятельства возникают как при попытках обнаружения и определения структуры привнесенных частиц, так и решении проблем синтеза слоистых структур для получения требуемых спектральных свойств (см. [99], [100]). Кроме того, детальный анализ свойств 3D резонаторов плазмонного нанолазера также требует полного учета взаимодействия слоистой частицы с активной средой и прозрачной подложкой (см. [17]). В этом плане следует отметить, что в современной литературе, кроме работ авторов [101], [102], отсутствуют результаты моделирования поведения полей слоистых частиц, расположенных на подложке, с учетом эффекта нелокальности.

Перейдем к построению математической модели. Пусть все пространство ${{\mathbb{R}}^{3}}$ разделено на два полупространства: активная среда ${{D}_{0}}:(z > 0)$ и прозрачная подложка ${{D}_{1}}:(z < 0)$. Обозначим через ${{\Sigma }}:\left( {z = 0} \right)$ плоскую границу раздела сред. Пусть сферическая слоистая частица целиком располагается в верхнем полупространстве ${{D}_{0}}$. Внутреннее диэлектрическое ядро частицы обозначим ${{D}_{i}}$, а область металлической оболочки ${{D}_{s}}$. Соответствующие сферические поверхности будем обозначать как $\partial {{D}_{{i,s}}}$. Все среды предполагаются немагнитными, а их диэлектрические проницаемости обозначим ${{\varepsilon }_{\nu }},\;\nu = 0,~1,~i,~s$.

Перейдем к формулировке математической постановки граничной задачи рассеяния для системы Максвелла с учетом модели ОНО. Она во многом сходна с граничной задачей (4), но имеет существенные отличительные особенности. Обозначим $\{ {{{\mathbf{E}}}^{0}},{{{\mathbf{H}}}^{0}}\} $ – поле плоской электромагнитной волны линейной поляризации, распространяющейся из верхнего полупространства в полуплоскости $\varphi = \pi $ под углом $\pi - {{\theta }_{0}}$ относительно нормали к подложке, совпадающей с осью 0z. Тогда постановка граничной задачи рассеяния с учетом ЭН может быть записана в следующем виде:

с условиями излучения на бесконечности вида Здесь $\{ {{{\mathbf{E}}}_{\nu }},{{{\mathbf{H}}}_{\nu }}\} $ – полные поля в областях ${{D}_{\nu }},~\;\nu = 0,~1,~i$, $\{ {\mathbf{E}}_{\nu }^{s},{\mathbf{H}}_{\nu }^{s}\} $ – рассеянное поле в ${{D}_{{0,1}}}$, ${{{\mathbf{n}}}_{{i,s}}}$ – единичные внешние нормали к поверхностям $\partial {{D}_{{i,s}}}$, ${{\varepsilon }_{\nu }},\;\nu = 0,~1,~i,s$ – диэлектрические проницаемости сред в соответствующих областях, при этом $\operatorname{Im} {{\varepsilon }_{{i,0,1}}}$ = 0, $\operatorname{Im} {{\varepsilon }_{{s,L}}} \leqslant 0$, $k = \frac{\omega }{c}$, временная зависимость величин, как и прежде, выбрана в виде $\exp (j\omega t)$.

Отметим, что условия излучения выбраны таким образом, чтобы обеспечить стремление к нулю потока энергии на бесконечности в отсутствиe внешнего возбуждения (см. [103]).

Конкретизируем остальные величины, входящие в постановку задачи (5). Поле $\{ {\mathbf{E}}_{0}^{0},{\mathbf{H}}_{0}^{0}\} $ представляет собой результат решения задачи отражения и преломления поля плоской волны $\{ {{{\mathbf{E}}}^{0}},{{{\mathbf{H}}}^{0}}\} $ на плоскости раздела полупространств Σ. Рассеянные поля $\{ {\mathbf{E}}_{\nu }^{s},{\mathbf{H}}_{\nu }^{s}\} $ определяются как ${\mathbf{E}}_{\nu }^{s} = {{{\mathbf{E}}}_{\nu }} - {\mathbf{E}}_{\nu }^{0}$, ${\mathbf{H}}_{\nu }^{s} = {{{\mathbf{H}}}_{\nu }} - {\mathbf{H}}_{\nu }^{0}$, $\nu = 0,~1$. В силу построения поля внешнего возбуждения и граничных условий на Σ рассеянное поле $\{ {\mathbf{E}}_{{0,1}}^{s},{\mathbf{H}}_{{0,1}}^{s}\} $ также должно удовлетворять условиям сопряжения для тангенциальных компонент на бесконечной границе Σ.

Так как частица располагается целиком в верхнем полупространстве ${{D}_{{0,}}}$, то суммарное поле падающей и отраженной волн $\{ {\mathbf{E}}_{0}^{0},{\mathbf{H}}_{0}^{0}\} $ для P и S поляризаций приобретает вид

Здесь ${\mathbf{E}}_{0}^{{P,S\left( - \right)}} = {{{\mathbf{E}}}^{{0\left( {P,S} \right)}}}$, ${\mathbf{H}}_{0}^{{P,S\left( - \right)}} = {{{\mathbf{H}}}^{{0\left( {P,S} \right)}}}$, ${\mathbf{E}}_{0}^{{P\left( \pm \right)}} = {\mathbf{e}}_{0}^{ \pm } \cdot \chi _{0}^{ \pm }~$, ${\mathbf{H}}_{0}^{{P\left( \pm \right)}} = - {{{\mathbf{e}}}_{y}}{{n}_{0}} \cdot \chi _{0}^{ \pm }$, ${\mathbf{E}}_{0}^{{S\left( \pm \right)}} = {{{\mathbf{e}}}_{y}} \cdot \chi _{0}^{ \pm }$, ${\mathbf{H}}_{0}^{{S\left( \pm \right)}} = {\mathbf{e}}_{0}^{ \pm }{{n}_{0}} \cdot \chi _{0}^{ \pm }$, $\chi _{0}^{ \pm } = \exp \{ - j{{k}_{0}}(x\sin {{\theta }_{0}} \pm z\cos {{\theta }_{0}})\} $, ${{n}_{0}} = \sqrt {{{\varepsilon }_{0}}} $, ${{k}_{0}} = k{{n}_{0}}$, ${\mathbf{e}}_{0}^{ \pm } = ( \mp {{{\mathbf{e}}}_{x}}\cos {{\theta }_{0}} + {{{\mathbf{e}}}_{z}}\sin {{\theta }_{0}})$.

Поскольку мы будем рассматривать наряду с возбуждением плоской волной, падающей из ${{D}_{0}}$, также случай возбуждения неизлучающей волной, распространяющейся из призмы, то в этом случае преломленное поле в верхнее полупространство имеют вид

(в данном случае в отличиe от падения сверху), ${{\theta }_{0}}$ – угол преломления, под которым волна проходит из призмы в ${{D}_{0}}$, а ${{R}_{{P,S}}};~\;{{T}_{{P,S}}}~$ – коэффициенты отражения и преломления Френеля (см. [104]).

Построим приближенное решение задачи (4) для рассеянного поля в ${{D}_{0}}$ с учетом осевой симметрии и поляризации, удовлетворяя квазиклассической системе уравнений Максвелла во всех областях постоянства параметров среды, условиям излучения и дополнительно условиям сопряжения для полей на бесконечной границе раздела полупространств $\Sigma $. В основу представления для рассеянного поля положим компоненты Фурье тензора Грина полупространства, которые могут быть записаны в виде интегральных представлений Вейля–Зоммерфельда (см. [50])

где ${{J}_{m}}({\kern 1pt} .{\kern 1pt} )$ – цилиндрическая функция Бесселя, точка $\zeta = (\rho ,z)$ располагается в полуплоскости $\varphi $ = const, а точки локализации мультиполей распределены вдоль оси симметрии ${{z}_{n}} \in 0z$ строго внутри ${{D}_{i}} \cup {{D}_{s}}$. Спектральные функции электрического и магнитного типов ${v}_{{11}}^{{e,h}},\;{v}_{{31}}^{{e,h}}$ обеспечивают выполнение условий сопряжения на границе интерфейса $z = 0$. В данном случае для них справедливы следующие выражения: Спектральные коэффициенты A, B определяются из одномерной задачи с условиями сопряжения при $z = 0$, откуда легко получается, что где ${{\eta }_{\nu }} = \sqrt {{{\lambda }^{2}} - k_{\nu }^{2}} $, $\gamma _{\nu }^{e} = {{\eta }_{\nu }}$, $\gamma _{\nu }^{h} = {{\eta }_{\nu }}{\text{/}}{{\varepsilon }_{\nu }}$, $\nu = 0,1.$

Будем строить решения с учетом поляризации. Начнем с P поляризации. Тогда для построения приближенного решения для рассеянного поля в ${{D}_{0}}$ мы будем использовать векторные потенциалы, которые в цилиндрической системе координат могут быть записаны как

Само же приближенное решение принимает вид Последний член в формуле (5.8) соответствует вертикальным электрическим диполям. Отметим, что представление для решения в областях ${{D}_{{0,1}}}$ записывается с единым набором амплитуд $\{ p_{{mn}}^{{0,P}},q_{{mn}}^{{0,P}},r_{n}^{{0,P}}\} $. Последнее обстоятельство является следствием использования тензора Грина для представления решения (5.8).

Обратимся к случаю S поляризации. В этом случае соответствующие векторные потенциалы принимают вид

Соответственно приближенное решение для S поляризации будет В данном случае последнее слагаемое (5.9) соответствует вертикальным магнитным диполям. Итак, приближенные решения (5.8), (5.9) удовлетворяют системе уравнений Максвелла (5.1) в ${{D}_{{0,1}}}$, условиям излучения, а также условиям сопряжения (5.5) на бесконечной границе раздела полупространств в силу использования компонент тензора Грина. Таким образом, для завершения построения приближенного решения остается построить приближенное решение для внутренней области. Заметим, что это построение нами уже представлено в предыдущем разделе (4.11), (4.13).

Подчеркнем, что численный алгоритм определения неизвестных амплитуд ДИ остается в точности тем же самым, что и в случае свободного пространства. Таким образом, используются проекционные соотношения вида (4.15), (4.16).

Для вычисления характеристик рассеяния в дальней зоне нам понадобится диаграмма направленности рассеянного поля ${{{\mathbf{F}}}^{{0,1}}}(\theta ,\varphi )$, которая определяется в верхнем и нижнем полупространствах ${{D}_{{0,1}}}$ как

Тогда компоненты $(\theta ,\varphi )$ диаграммы на ${{{{\Omega }}}^{ + }} = \{ 0^\circ \leqslant \theta \leqslant 90^\circ ;\;0^\circ \leqslant \varphi \leqslant 360^\circ \} $ – единичной верхней полусфере для Р поляризации – принимают вид где соответствующие спектральные функции $\bar {G}_{n}^{{e,h(0)}},\;\bar {g}_{n}^{{e,h(0)}}$ могут быть представлены как Отметим, что при ${{\varepsilon }_{0}} = {{\varepsilon }_{1}}$ представления для диаграмм рассеяния (5.11), (5.12) переходят в соответствующие представления для свободного пространства (4.18), (4.19).

Используя технику, аналогичную [50], нетрудно получить представления для диаграммы направленности в нижнем полупространстве ${{{{\Omega }}}^{ - }} = \{ 90^\circ \leqslant \theta \leqslant 180^\circ ;\;0^\circ \leqslant \varphi \leqslant 360^\circ \} $, компоненты которой приобретают вид

и, соответственно, для S поляризации где спектральные функции $\bar {G}_{n}^{{e,h(1)}},\;\bar {g}_{n}^{{e,h(1)}}$ записываются в виде Таким образом, определяя амплитуды ДИ для рассеянного поля, из проекционных соотношений (4.15), (4.16) можно легко вычислять компоненты диаграммы направленности (5.11)–(5.14) всюду на единичной сфере, а также поля (5.8), (5.9) в непосредственной близости от оболочки слоистой частицы. Подчеркнем еще раз, что диаграмма рассеяния, как в нижнем, так и верхнем полупространстве, вычисляется на основе одних и тех же амплитуд ДИ $\{ p_{{mn}}^{0},~q_{{mn}}^{0},r_{n}^{0}\} $, что является следствием использования тензора Грина, реализующего единое представление для рассеянного поля всюду в ${{D}_{{0,1}}}$.

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

Будем рассматривать слоистую частицу, состоящую из SiO₂ (n_i = 1.46) ядра сферической формы, покрытую золотой (Au) оболочкой толщиной d, располагающейся в активной среде R6G (n_e = 1.326) на поверхности стеклянной призмы ВК7 (n₁= 1.52). Нас интересуют численные результаты анализа коэффициента усиления ближнего поля (4.23) и сечения поглощения на внешней границе слоя в диапазоне длин волн

Поскольку мы будем рассматривать возбуждение плоской волной, распространяющейся не только из верхнего полупространства, но и из-под поверхности призмы, то необходимо определить область неизлучающих волн, располагающуюся за углом полного внутреннего отражения. В соответствии с законом Снеллиуса (см. [104]) имеет место соотношение ${{n}_{0}}\sin {{\theta }_{0}} = {{n}_{1}}\sin {{\theta }_{1}}$. Откуда преломленный угол ${{\theta }_{0}} = \arcsin [({{n}_{1}}{\text{/}}{{n}_{0}})\sin {{\theta }_{1}}]$. В случае, когда волна падает из более плотной среды в менее плотную ${{n}_{1}} > {{n}_{0}}$, существует критический угол ${{\theta }_{c}} = \arcsin ({{n}_{0}}{\text{/}}{{n}_{1}}),$ за которым $({{\theta }_{1}} > {{\theta }_{c}})$ волна не проходит в верхнее полупространство, так как полностью отражается от поверхности призмы $\Sigma $. При этом энергия распространяется вдоль поверхности раздела полупространств и экспоненциально затухает в направлении, перпендикулярном к поверхности призмы. В этом случае оказывается, что $\sin {{\theta }_{0}} > 1$ и $\cos {{\theta }_{0}} = - j\sqrt {{{{\sin }}^{2}}{{\theta }_{0}} - 1} $. Таким образом, амплитуда плоской волны в ${{D}_{0}}$ приобретает вид $\exp \{ - j{{k}_{0}}x\sin {{\theta }_{0}}\} \exp \left\{ { - {{k}_{0}}z{\text{ }}\sqrt {{{{\sin }}^{2}}{{\theta }_{0}} - 1} } \right\}$. В нашем конкретном случае критический угол равняется ${{\theta }_{c}} = 60.735^\circ .$ Следовательно, при углах падения ${{\theta }_{1}} > 60.735^\circ $ мы попадаем в область неизлучающих волн.

Зафиксируем диаметр диэлектрического ядра слоистой частицы ядра ${{D}_{c}} = 14$ нм, а толщины золотой оболочки будут, как и ранее, $d = 2.5,~\;3.5$ нм. Квантовые параметры будем выбирать в соответствии с (4.24).

На фиг. 5а приведены значения коэффициента усиления (4.23), соответствующего возбуждению плоской волной, распространяющейся из верхнего полупространства под углом ${{\theta }_{0}} = 45^\circ $. Как показали проведенные исследования (см. [101]), этот угол реализует максимальные значения интенсивности ближнего поля. На фиг. 5а видно, что уменьшение толщины золотой пленки приводит к смещению ПР в область длинных волн, при этом учет ЭН приводит к снижению амплитуды ПР до 45%. Сходные с предыдущим результаты можно видеть на фиг. 5б, где изображено сечение поглощения. Однако смещение ПР не сопровождается в данном случае возрастанием его амплитуды, кроме того, учет ЭН приводит к менее заметному уменьшению амплитуды ПР.

На фиг. 6а приведены результаты возбуждения частицы неизлучающей волной, соответствующей углу падения ${{\theta }_{1}} = 62^\circ $. Несмотря на сходство с фиг. 5а в этом случае амплитуда ПР существенно больше. Что касается снижения при учете ЭН, то оно близко к отмеченному на фиг. 5а. Последняя кривая соответствует толщине слоя $d = 1.5$ нм с учетом ЭН и нужна для интерпретации результатов, относящихся к асимметричным частицам. На фиг. 6б изображено сечение поглощения и результаты качественно схожи с изображенными на фиг. 5б. Однако легко заметить, что величины сечения поглощения $\sigma _{{abs}}^{P}$ возросли почти в 5 раз по сравнению с возбуждением частицы, распространяющейся волной.

Далее мы будем исследовать влияния асимметрии на оптические характеристики. На фиг. 7 приведены значения коэффициента усиления при возбуждении распространяющейся волной для толщины пленки $d = 2.5$ нм при различных величинах смещения sh центра ядра относительно центра внешней поверхности оболочки по оси симметрии. Значение смещения ± соответствуют смещению центра вверх и вниз соответственно. Видно, что при любом смещении ПР сдвигается вправо и возрастает. Это еще более заметно на фиг. 8 с результатами, соответствующими учету ЭН: любое смещение вверх и вниз ведет к смещению вправо и увеличению амплитуды ПР.

Сходные с фиг. 7 результаты имеют место для случая возбуждения частицы неизлучающей волной (фиг. 9). Снова любая асимметрия ведет к увеличению амплитуды ПР, что легко объяснить, рассматривая результат на фиг. 6а для $d = 1.5$ нм. Однако в случае сечения поглощения наблюдается противоположная тенденция, что можно видеть на фиг. 10, где амплитуда ПР снижается вместе с увеличением асимметрии.

Резюмируя полученные результаты, отмечаем, что асимметрия слоистой частицы приводит к изменениям в положении ПР и его амплитуде. При этом в случае с коэффициентом усиления асимметрия ведет к его увеличению, в то время как в случае с сечением поглощения она ведет к снижению величины. По-видимому, увеличение коэффициента усиления связано с появлением тонкой части оболочки, что, как отмечено на фиг. 6а, приводит к существенному росту коэффициента усиления. Во всех рассмотренных случаях учет эффекта нелокальности приводит к снижению ПР при небольшом сдвиге влево в области длин волн (blue shift).

ЗАКЛЮЧЕНИЕ

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