1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 63, № 11, 2023

Журнал вычислительной математики и математической физики, 2023, T. 63, № 11, стр. 1911-1921

Анализ влияния квантовых эффектов на оптические характеристики плазмонных наночастиц методом дискретных источников

Ю. А. Еремин 1, *, В. В. Лопушенко 1, **

1 МГУ им. М.В. Ломоносова, ВМК
119991 Москва, Ленинские горы, 1, стр. 52, Россия

* E-mail: eremin@cs.msu.ru
** E-mail: lopushnk@cs.msu.ru

Поступила в редакцию 24.05.2023
После доработки 24.05.2023
Принята к публикации 25.07.2023

Полный текст (PDF)

Аннотация

Метод дискретных источников адаптирован для исследования проявления поверхностных квантовых эффектов, основанных на мезоскопических граничных условиях с параметрами Фейбельмана. Проводится сравнительный анализ влияния объемных эффектов нелокальности в рамках теории обобщенного нелокального отклика и поверхностных эффектов на оптические характеристики золотых и серебряных наночастиц. Установлено, что если учет эффекта нелокальности для благородных металлов всегда сопровождается снижением амплитуды и сдвигом в коротковолновую область, то влияние поверхностного эффекта существенно зависит от геометрии частиц. При этом мезоскопические граничные условия в значительной степени восстанавливают амплитуду плазмонного резонанса по сравнению с объемным эффектом нелокальности. Это различие особенно заметно при сравнительном анализе коэффициента усиления поля на поверхности частиц. Также установлено существенное отличие в поведении плазмонного резонанса для золотых и серебряных частиц для случая мезоскопических граничных условий. Библ. 27. Фиг. 6.

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

1. ВВЕДЕНИЕ

Классическая теория Максвелла применительно к исследованию оптических свойств металлических структур нашла широкое применение в плазмонике, способствуя успешному описанию особенностей взаимодействия света с веществом в масштабах выше 10 нм [1]. Углубленное исследование наноплазмонных структур позволило обнаружить множество оптических эффектов, одним из которых является локализованный поверхностный плазмонный резонанс (ПР). ПР способен на порядки усиливать поля и концентрировать энергию в объемах, намного превосходящих рэлеевский оптический предел. Подобные эффекты имеют фундаментальное значение для развития новых технологий, включая обнаружение и терапию онкологических образований, исследования свойств отдельных молекул, разработку резонаторов плазмонного нанолазера, а также создание конгломератов, позволяющих концентрировать и удерживать солнечную энергию [2].

Вместе с тем современное развитие прорывных технологий нанопроизводства ведет к необходимости манипулировать наноразмерными структурами размерами 1–10 нм [3]. В таком масштабе описание происходящих процессов посредством классической теории Максвелла оказывается недостаточным, так как она дает лишь качественную картину наблюдаемых явлений. Это связано прежде всего с возникновением таких квантовых эффектов, как нелокальность и концентрация свободных электронов вблизи и внутри плазмонного металла [4, 5]. Следует отметить, что существующие чисто квантовомеханические подходы практически неприменимы для полного описания подобных эффектов в плазмонном материале реальных наноструктур [5]. Например, даже после значительных приближений и упрощений очень хорошо зарекомендовавшая себя теория функционала плотности, зависящей от времени (TDDFT) [6], не позволяет исследовать металлические наночастицы размером более нескольких нанометров [7].

В связи с этим было предложено множество методов, основанных на квазиклассических теориях, для описания квантовых эффектов в целях исследования плазмонных объектов мезоскопического масштаба с характерной длиной от 1 до 10 нм [8]. Среди этих теорий наибольшее распространение получила гидродинамическая теория Друде [4] и ее обобщение – теория обобщенного нелокального отклика (GNOR) [9]. В частности, GNOR вместе с граничными условиями непротекания для тока на границе раздела металл–диэлектрик способна учитывать эффект нелокальности и поверхностное затухание Ландау. Однако даже теория GNOR, учитывающая эффект нелокальности, не описывает выброс электронов за пределы металла [10].

В последнее время был разработан подход, позволяющий учитывать квантовые эффекты в диапазоне 1–10 нм, так называемом мезоскопическом масштабе. Он оперирует функциями поверхностного отклика (SRF) – d-параметрами Фейбельмана ${{d}_{ \bot }}$ и ${{d}_{\parallel }}$ [11, 12]. При этом вещественные части ${{d}_{ \bot }}$ и ${{d}_{\parallel }}$ описывают положение центра масс индуцированного заряда и плотности тока вблизи границы металла соответственно, тогда как их мнимые части вносят вклад в поверхностное поглощение и затухание. Использование функций поверхностного отклика направлено на то, чтобы учесть доминирующие квантовые явления вблизи поверхности, при этом рассмотрение взаимодействия света и вещества в объеме плазмонного металла производится в рамках классической теории Максвелла [13, 14]. Недавно $d$-параметры были успешно включены в граничные условия для электромагнитных полей с целью учета доминирующего неклассического поверхностного оптического отклика [13]. Этот подход имеет большое значение при изучении плазмонных систем мезоскопического масштаба, поскольку он значительно повышает эффективность вычислений по сравнению с трудоемким подходом TDDFT [15, 16]. Следует отметить, что функции поверхностного отклика представляются одним из многообещающих путей преодоления разрыва между чисто квантово-механическими подходами и феноменологическими моделями квантовой наноплазмоники [17, 18].

2. МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА ЗАДАЧИ В РАМКАХ ТЕОРИИ ОБОБЩЕННОГО НЕЛОКАЛЬНОГО ОПТИЧЕСКОГО ОТКЛИКА (GNOR)

Будем рассматривать задачу дифракции поля плоской электромагнитной волны $\left\{ {{{{\mathbf{E}}}_{0}},{{{\mathbf{H}}}_{0}}} \right\}$ на однородной осесимметричной плазмонной наночастице с внутренней областью ${{D}_{i}}$ и внешней ${{D}_{e}}$, разделенных гладкой замкнутой поверхностью $\partial {{D}_{i}}$. Будем предполагать, что все среды немагнитные, а зависимость от времени выбрана в виде ${\text{exp}}\left( {j\omega t} \right)$.

Суть теории нелокальности состоит в том, что локальная связь между электрическим полем и смещением ${{{\mathbf{D}}}_{i}}\left( {{\mathbf{r}},{{\omega }}} \right) = {{{{\varepsilon }}}_{i}}\left( {{\omega }} \right){{{\mathbf{E}}}_{i}}\left( {{\mathbf{r}},{{\omega }}} \right)$ преобразуется в интегральную ${{{\mathbf{D}}}_{i}}\left( {{\mathbf{r}},{{\omega }}} \right) = \int {{{{{\varepsilon }}}_{i}}\left( {{\mathbf{r}},{\mathbf{r}}{\kern 1pt} ',{{\omega }}} \right){{E}_{i}}\left( {{\mathbf{r}}{\kern 1pt} ',{{\omega }}} \right){{d}^{3}}{\mathbf{r}}{\kern 1pt} '} $. Последнее обстоятельство влечет за собой возникновение в ${{D}_{i}}$ продольного поля ${\mathbf{E}}_{i}^{L}\,:\,{\text{\;rot}}{\mathbf{E}}_{i}^{L} = 0$ вдобавок к классическому поперечному ${\mathbf{E}}_{i}^{T}\,:\,{\text{div}}{\mathbf{E}}_{i}^{T} = 0$. Таким образом, поле внутри частицы ${{D}_{i}}$ представляется в виде суммы поперечного и продольного полей ${{{\mathbf{E}}}_{i}} = {\mathbf{E}}_{i}^{T} + {\mathbf{E}}_{i}^{L}$. Можно показать, что эти поля внутри ${{D}_{i}}$ удовлетворяют следующим уравнениям Гельмгольца:

(1а)
$\Delta {\mathbf{E}}_{i}^{T}\left( {{\mathbf{r}},{{\omega }}} \right) + k_{T}^{2}\left( {{\omega }} \right){\mathbf{E}}_{i}^{T}\left( {{\mathbf{r}},{{\omega }}} \right) = 0,\,\,\,\,\Delta {\mathbf{E}}_{i}^{L}\left( {{\mathbf{r}},{{\omega }}} \right) + k_{L}^{2}\left( {{\omega }} \right){\mathbf{E}}_{i}^{L}\left( {{\mathbf{r}},{{\omega }}} \right) = 0.$

Во внешней области рассеянное поле $\left\{ {{{{\mathbf{E}}}_{e}},{{{\mathbf{H}}}_{e}}} \right\}$ удовлетворяет классической системе уравнений Максвелла

(1б)
${\text{rot}}{{{\mathbf{H}}}_{e}} = jk{{\varepsilon }_{e}}{{{\mathbf{E}}}_{e}};\,\,\,\,{\text{\;rot}}{{{\mathbf{E}}}_{e}} = - jk{{{\mathbf{H}}}_{e}};\,\,{\text{в}}\,\,{{D}_{e}}.$

На поверхности частицы имеют место условия сопряжения для тангенциальных компонент полей

(1в)
${{{\mathbf{n}}}_{i}} \times \left( {{{{\mathbf{E}}}_{i}}\left( P \right) - {{{\mathbf{E}}}_{e}}\left( P \right)} \right) = {{{\mathbf{n}}}_{i}} \times {{{\mathbf{E}}}_{0}}\left( P \right);P \in \partial {{D}_{i}},\,\,\,\,{{{\mathbf{n}}}_{i}} \times \left( {{{{\mathbf{H}}}_{i}}\left( P \right) - {{{\mathbf{H}}}_{e}}\left( P \right)} \right) = {{{\mathbf{n}}}_{i}} \times {{{\mathbf{H}}}_{0}}\left( P \right);$
с дополнительным граничным условием на поверхности металл–диэлектрик, которое формулируется как обращение в нуль нормальной компоненты тока проводимости и может быть преобразовано к виду [9]

(1г)
${{\varepsilon }_{b}}{{{\mathbf{n}}}_{i}} \cdot {{{\mathbf{E}}}_{i}}\left( Q \right) = {{\varepsilon }_{e}}{{{\mathbf{n}}}_{i}} \cdot {{{\mathbf{E}}}_{e}}\left( Q \right),\,\,\,\,Q \in \partial {{D}_{i}}.$

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

(1д)
$\mathop {\lim }\limits_{r \to \infty } r\left( {\sqrt {{{\varepsilon }_{e}}} {{{\mathbf{E}}}_{e}} \times \frac{{\mathbf{r}}}{r} - {{{\mathbf{H}}}_{e}}} \right) = 0;\,\,\,\,r = \left| M \right|.{\text{\;}}$

Здесь ${{\varepsilon }_{{i,e}}}$ – диэлектрические проницаемости сред в соответствующих областях, при этом ${\text{Im}}{{\varepsilon }_{e}}$ = 0, ${\text{Im}}{{\varepsilon }_{i}} \leqslant 0$, $k = \frac{\omega }{c}$, $k_{T}^{2}\left( \omega \right) = {{k}^{2}}{{\varepsilon }_{i}}\left( \omega \right)~$ – поперечное волновое число, а $k_{L}^{2}\left( \omega \right) = \frac{{{{\varepsilon }_{i}}\left( \omega \right)}}{{{{\xi }^{2}}\left( \omega \right)}}$ продольное. Соответствующие квантовые параметры имеют вид: $\xi ~$ – длина корреляции в рамках теории GNOR, ${{\xi }^{2}}\left( \omega \right) = {{\varepsilon }_{b}}\frac{{{{\beta }^{2}} + D\left( {\gamma + j\omega } \right)}}{{{{\omega }^{2}} - j\gamma \omega }}$, ${{\varepsilon }_{b}} = {{\varepsilon }_{i}} + \frac{{\omega _{p}^{2}}}{{{{\omega }^{2}} - j\gamma \omega }}$, ${{\omega }_{p}}$ – плазменная частота металла, ${{{{\beta }}}^{2}} = \left( {\frac{3}{5}} \right)v_{{\text{F}}}^{2}$, ${{v}_{{\text{F}}}}$ – скорость Ферма, ${{\gamma \;}}$ – скорость затухания Друде, $D{\text{\;}}$ – коэффициент диффузии электронов в металле. В случае, когда $D = 0$ теория GNOR переходит в гидродинамическую теорию Друде [4].

Таким образом, математическая постановка граничной задачи дифракции в рамках теории GNOR включает в себя уравнения Гельмгольца внутри ${{D}_{i}}$ (1а), систему уравнений Максвелла вне рассеивателя (1б), условия сопряжения для тангенциальных компонент полей на границе $\partial {{D}_{i}}$ (1в), плюс дополнительное граничное условие (1г) и условия излучения на бесконечности (1д). Мы будем полагать, что граничная задача (1) разрешима единственным образом. Для случая сферы это может быть установлено аналитически [4].

3. ГРАНИЧНАЯ ЗАДАЧА ДИФРАКЦИИ С МЕЗОСКОПИЧЕСКИМИ ГРАНИЧНЫМИ УСЛОВИЯМИ

В данном случае математическая постановка задачи дифракции может быть записана в следующем виде:

(2а)
${\text{rot}}{{{\mathbf{H}}}_{{i,e}}} = jk{{\varepsilon }_{{i,e}}}{{{\mathbf{E}}}_{{i,e}}};{\text{\;rot}}{{{\mathbf{E}}}_{{i,e}}} = - jk{{{\mathbf{H}}}_{{i,e}}}\,\,\,\,{\text{в}}\,\,{{D}_{{e,i}}};$
(2б)
$\begin{gathered} {{{\mathbf{n}}}_{i}} \times \left( {{{{\mathbf{E}}}_{i}}\left( P \right) - {{{\mathbf{E}}}_{e}}\left( P \right) - {{{\mathbf{E}}}_{0}}\left( P \right)} \right) = - {{d}_{ \bot }}{{{\mathbf{n}}}_{i}} \times \nabla {\text{\;}}\left\{ {{{{\mathbf{n}}}_{i}} \cdot \left( {{{{\mathbf{E}}}_{i}}\left( P \right) - {{{\mathbf{E}}}_{e}}\left( P \right) - {{{\mathbf{E}}}_{0}}\left( P \right)} \right)} \right\}, \\ {{{\mathbf{n}}}_{i}} \times \left( {{{{\mathbf{H}}}_{i}}\left( P \right) - {{{\mathbf{H}}}_{e}}\left( P \right) - {{{\mathbf{H}}}_{0}}\left( P \right)} \right) = - j{{\omega }}{{d}_{\parallel }}\left\{ {{{{\mathbf{n}}}_{i}} \times \left[ {{{{\mathbf{D}}}_{i}}\left( P \right) - {{{\mathbf{D}}}_{e}}\left( P \right) - {{{\mathbf{D}}}_{0}}\left( P \right)} \right]} \right\} \times {{{\mathbf{n}}}_{i}},~ \\ P \in \partial {{D}_{i}}; \\ \end{gathered} $
(2в)
$\mathop {\lim }\limits_{r \to \infty } r\left( {\sqrt {{{\varepsilon }_{e}}} {{{\mathbf{E}}}_{{\text{e}}}} \times \frac{{\mathbf{r}}}{r} - {{{\mathbf{H}}}_{e}}} \right) = 0;\,\,\,\,r = \left| M \right|.{\text{\;}}$

Здесь ${{d}_{ \bot }},{\text{\;}}~{{d}_{\parallel }}$ – параметры Фейбельмана (ПФ). Формально они определяются следующим образом:

${{d}_{ \bot }}\left( {{\omega }} \right) = \frac{{\int {\rho \left( {x,\omega } \right)xdx} }}{{\int {\rho \left( {x,\omega } \right)dx} }},{\text{\;}}\,\,\,{{d}_{\parallel }}\left( {{\omega }} \right) = {\text{\;}}\frac{{\int {{{\partial }_{x}}{{J}_{y}}\left( {x,\omega } \right)xdx} }}{{\int {{{\partial }_{x}}{{J}_{y}}\left( {x,\omega } \right)dx} }}{\text{\;}},$
где $\rho \left( {x,\omega } \right)$ – плотность индуцированных поверхностных зарядов, а ${{J}_{y}}\left( {x,\omega } \right)$ – токов. Уместно отметить, что для случая границы раздела сред металл–диэлектрик параметр Фейбельмана ${{d}_{\parallel }} = 0$ [15]. Параметр ${{d}_{ \bot }}$ представляет собой комплексную функцию, зависящую от длины падающей волны. Его действительная часть описывает положение центроида индуцированных зарядов по отношению к границе, а мнимая отвечает за поглощение энергии вблизи поверхности (затухание Ландау) [18].

Как и в предыдущем случае, будем полагать, что граничная задача (2) имеет единственное решение. Для случая однородной сферы это установлено в [14].

4. МЕТОД ДИСКРЕТНЫХ ИСТОЧНИКОВ

Для построения приближенных решений граничных задач (1) и (2) мы используем метод дискретных источников (МДИ) [19, 20], который представляет собой строгий полуаналитический поверхностно ориентированный метод. В рамках МДИ представление для поля строится как конечная линейная комбинация распределенных мультиполей низшего порядка [21], удовлетворяющая уравнениям Максвелла (Гельмгольца) везде вне разрывов среды и условиям излучения Сильвера-Мюллера на бесконечности в явном аналитическом виде. Соответствующие амплитуды дискретных источников определяются из граничных условий, заданных на поверхности частицы. Основное преимущество МДИ состоит в том, что он позволяет оценивать погрешность полученного решения путем вычисления невязки полей на поверхности частицы. Последнее обстоятельство позволяет с заданной точностью рассчитывать оптические характеристики полей вблизи поверхности частиц. Кроме того, мы хотели бы подчеркнуть, что МДИ весьма удобен для анализа задачи рассеяния с мезоскопическими граничными условиями (2), потому что не требуется никаких дополнительных усилий для учета дифференциального оператора в граничных условиях. Так как в рамках МДИ поля вблизи поверхности частицы представляют собой аналитические функции, можно вычислять производные любого типа и порядка.

Начнем с построения приближенного решения граничной задачи (1) для случая р‑поляризации падающей плоской волны. Именно в этом случае возникает наиболее заметный ПР [20]. Поле внешнего возбуждения может быть записано в виде

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

Будем строить поперечные поля в областях ${{D}_{{i,e}}}$ на основе векторных потенциалов, индуцированных источниками, распределенными вдоль оси симметрии рассеивателя. Потенциалы могут быть записаны, см. [22], в виде

(4)
$\begin{gathered} {\mathbf{A}}_{{mn}}^{{1,\alpha }} = Y_{m}^{\alpha }\left( {{{\zeta }},z_{n}^{\alpha }} \right)\left\{ {{{{\mathbf{e}}}_{\rho }}{\text{cos}}\left[ {\left( {m + 1} \right){{\varphi }}} \right] - {{{\mathbf{e}}}_{{{\varphi }}}}{\text{sin}}\left[ {\left( {m + 1} \right){{\varphi }}} \right]} \right\};\,\,\,\,{{\alpha }} = i,e; \\ {\mathbf{A}}_{{mn}}^{{2,\alpha }} = Y_{m}^{\alpha }\left( {{{\zeta }},z_{n}^{\alpha }} \right)\left\{ {{{{\mathbf{e}}}_{\rho }}{\text{sin}}\left[ {\left( {m + 1} \right)\varphi } \right] + {{{\mathbf{e}}}_{{{\varphi }}}}{\text{cos}}\left[ {\left( {m + 1} \right){{\varphi }}} \right]} \right\};\,\,\,\,{\mathbf{A}}_{n}^{{3,{{\alpha }}}} = Y_{0}^{\alpha }\left( {{{\zeta }},z_{n}^{\alpha }} \right){{{\mathbf{e}}}_{z}}. \\ \end{gathered} $

Здесь соответствующие функции имеют вид

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

В свою очередь, приближенное решение для продольных полей внутри частицы ${{D}_{i}}{\text{\;}}$строится на основе скалярных потенциалов, которые удовлетворяют уравнению Гельмгольца вида (1а):

(5)
$\Psi _{{mn}}^{P} = j_{{m + 1}}^{{}}\left( {{{k}_{L}}{{R}_{{z_{n}^{i}}}}} \right){{\left( {\frac{\rho }{{{{R}_{{z_{n}^{i}}}}}}} \right)}^{{m + 1}}}\cos \left[ {\left( {m + 1} \right)\varphi } \right],\,\,\,\,~\Psi _{n}^{{}} = {{j}_{0}}\left( {{{k}_{L}}{{R}_{{z_{n}^{i}}}}} \right).$

Тогда конкретные представления для поперечных и продольных полей в случае р‑поляризации приобретают следующий вид:

(6)
$\begin{array}{*{20}{c}} {{\mathbf{E}}_{{T,\alpha }}^{N} = \mathop \sum \limits_{{\text{m}} = 0}^{\text{M}} \mathop \sum \limits_{n = 1}^{N_{\alpha }^{m}} \left\{ {p_{{mn}}^{\alpha }\frac{j}{{k{{\varepsilon }_{\alpha }}}}{\text{rot}}{\kern 1pt} {\text{rot}}{\kern 1pt} {\mathbf{A}}_{{mn}}^{{1,\alpha }} + q_{{mn}}^{\alpha }\frac{j}{{{{\varepsilon }_{\alpha }}}}{\text{rot}}{\kern 1pt} {\mathbf{A}}_{{mn}}^{{2,\alpha }}} \right\} + \mathop \sum \limits_{n = 1}^{N_{\alpha }^{0}} r_{n}^{\alpha }\frac{j}{{k{{\varepsilon }_{\alpha }}}}{\text{rot}}{\kern 1pt} {\text{rot}}{\kern 1pt} {\mathbf{A}}_{n}^{{3,\alpha }};} \\ {{\mathbf{E}}_{L}^{N} = \mathop \sum \limits_{{\text{m}} = 0}^{\text{M}} \mathop \sum \limits_{n = 1}^{N_{L}^{m}} p_{{mn}}^{L}\nabla \Psi _{{mn}}^{P} + \mathop \sum \limits_{n = 1}^{N_{L}^{0}} r_{n}^{L}\nabla {{\Psi }}_{{\text{n}}}^{{}}{\text{\;}};~\,\,\,\,{\mathbf{H}}_{\alpha }^{N} = \frac{j}{k}{\text{rot}}{\mathbf{E}}_{\alpha }^{{~N}}{\text{\;}},{\text{\;}}\,\,\,\alpha = i,e.~} \end{array}$

Легко убедиться, что построенные поля удовлетворяют всем условиям граничной задачи (1) за исключением граничных условий (1в), (1г). Неизвестные амплитуды дискретных источников ${\mathbf{p}}_{m}^{N} = \left\{ {p_{{mn}}^{{e,i}},q_{{mn}}^{{e,i}},r_{n}^{{e,i}};p_{{mn}}^{L},r_{n}^{L}} \right\}$ как раз и определяются из этих граничных условий.

Что касается приближенного решения граничной задачи (2), то оно имеет вид (6), но без продольного поля, т.е. остаются только представления для полей $\left\{ {{\mathbf{E}}_{{T,\alpha }}^{N},{\mathbf{H}}_{\alpha }^{N}} \right\},{\text{\;}}\alpha = e,i$.

Остановимся кратко на схеме численного алгоритма. Он включает в себя несколько этапов.

АЛГОРИТМ

Шаг 1. Проводится разложение в ряд Фурье по φ плоской волны $\chi \left( {{{\rho }},{{\varphi }},z} \right)$ и выбирается число членов ряда, достаточное для аппроксимации на поверхности.

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

Шаг 3. Для решения одномерной задачи аппроксимации используется обобщенный метод коллокаций [23]. Для этого выбираются точки на образующей, в которых сшиваются гармоники полей, причем число этих точек выбирается в 2–4 раза больше, чем суммарное число источников для каждой гармоники. Таким образом, матрицы оказываются переопределенными.

Шаг 4. Амплитуды определяются как нормальные псевдорешения соответствующих систем. При необходимости вычисляется невязка полей на поверхности в сеточной норме L2.

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

${{{{{\mathbf{E}}}_{e}}\left( M \right)} \mathord{\left/ {\vphantom {{{{{\mathbf{E}}}_{e}}\left( M \right)} {\left| {{{{\mathbf{E}}}^{0}}\left( M \right)} \right|}}} \right. \kern-0em} {\left| {{{{\mathbf{E}}}^{0}}\left( M \right)} \right|}} = \frac{{\exp \left( { - j{{k}_{e}}r} \right)}}{r}{\mathbf{F}}\left( {\theta ,\varphi } \right) + O({{r}^{{ - 2}}}),\,\,\,\,~r \to \infty .$

Следуя асимптотике рассеянного поля (6), компоненты $\theta ,\varphi $ диаграммы направленности для р-поляризации принимают следующий вид:

(7)
$\begin{gathered} F_{\theta }^{P}\left( {\theta ,\varphi } \right) = j{{k}_{e}}\mathop {{\kern 1pt} \sum }\limits_{m = 0}^M {{\left( {j{{k}_{e}}\sin \theta } \right)}^{m}}\cos \left( {m + 1} \right)\varphi {\kern 1pt} \mathop \sum \limits_{n = 1}^{N_{e}^{m}} \left\{ {p_{{mn}}^{e}\cos \theta + q_{{mn}}^{e}} \right\}\exp \left\{ {j{{k}_{e}}z_{n}^{e}\cos \theta } \right\} - \\ - \,\,j{{k}_{e}}\sin \theta {\kern 1pt} \mathop \sum \limits_{n = 1}^{N_{e}^{0}} r_{n}^{e}\exp \left\{ {j{{k}_{e}}z_{n}^{e}\cos \theta } \right\}, \\ F_{\varphi }^{P}\left( {\theta ,\varphi } \right) = - j{{k}_{e}}{\kern 1pt} \mathop \sum \limits_{m = 0}^M {{\left( {j{{k}_{e}}\sin \theta } \right)}^{m}}\sin \left( {m + 1} \right)\varphi {\kern 1pt} \mathop \sum \limits_{n = 1}^{N_{e}^{m}} \left\{ {p_{{mn}}^{e} + q_{{mn}}^{e}\cos \theta } \right\}\exp \left\{ {j{{k}_{e}}z_{n}^{e}\cos \theta } \right\}.~ \\ \end{gathered} $

Наш основной интерес будет связан с анализом поведения сечения экстинкции

(8)
$\sigma _{{{\text{ext}}}}^{P}\left( {{{\theta }_{0}}} \right) = - \frac{{4\pi }}{{{{k}_{e}}}}\operatorname{Im} F_{\theta }^{P}\left( {\pi - {{\theta }_{0}},\pi } \right),$
которое показывает, какую часть энергии плоской волны забирает частица на поглощение и рассеяние. Кроме того, мы будем вычислять коэффициент усиления поля на поверхности частицы $\partial {{D}_{i}}{\text{\;}},{\text{\;а\;именно}}$:

(9)
$F\left( \lambda \right) = \frac{{{{{\int\limits_{\partial {{D}_{i}}} {\left| {E_{e}^{N} + {{E}^{0}}} \right|} }}^{2}}d\sigma }}{{\int\limits_{\partial {{D}_{i}}} {{{{\left| {{{E}^{0}}} \right|}}^{2}}} d\sigma }}.~~$

Напомним, что МДИ позволяет вычислять ближние поля, компоненты диаграммы рассеяния и сечение экстинкции с учетом ЭН в аналитическом виде [25].

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

Перейдем к обсуждению результатов математического моделирования. Соответствующие квантовые параметры для золота и серебра в рамках модели GNOR будут выбираться равными [20]

Au: $\hbar {{\omega }_{p}} = 9.02\,\,{\text{eV}},\,\,\,\,\hbar \gamma = 0.071\,\,{\text{eV}},\,\,{{v}_{{\text{F}}}} = 1.39 \times {{10}^{{12}}}{{{{\mu м}}} \mathord{\left/ {\vphantom {{{{\mu м}}} {\text{c}}}} \right. \kern-0em} {\text{c}}},{\text{\;}}D = 8.62 \times {{10}^{8}}\,\,{{{{\mu }}{{{\text{м}}}^{2}}} \mathord{\left/ {\vphantom {{{{\mu }}{{{\text{м}}}^{2}}} {\text{c}}}} \right. \kern-0em} {\text{c}}}.$

Ag: $\hbar {{\omega }_{p}} = 8.99\,\,{\text{eV}},\,\,\,~\hbar \gamma = 0.025\,\,{\text{eV}},~\,\,\,\,{{v}_{{\text{F}}}} = 1.39 \times {{10}^{{12}}}\,\,\,\,{{\mu {\text{м}}} \mathord{\left/ {\vphantom {{\mu {\text{м}}} {\text{c}}}} \right. \kern-0em} {\text{c}}},\,\,\,\,~D = 9.62 \times {{10}^{8}}\,\,{{\mu {{{\text{м}}}^{2}}} \mathord{\left/ {\vphantom {{\mu {{{\text{м}}}^{2}}} {\text{c}}}} \right. \kern-0em} {\text{c}}}.$

Несмотря на то что теория GNOR адекватно описывает происходящие явления и ее результаты сопоставимы с измерениями и другими подходами [10], сами авторы теории считают, что она нуждается в некоторой коррекции [26]. В этой работе исследуется связь коэффициента диффузии с индуцированными поверхностными зарядами, описываемыми ПФ. В результате было получено следующее соотношение:

(10)
$D\left( \omega \right) = 4\sqrt 2 \beta ~\operatorname{Im} \left( {{{d}_{ \bot }}\left( \omega \right)} \right).$

Отсюда непосредственно следует, что коэффициент диффузии, который изначально принимал постоянное значение, может зависеть от длины волны. Мы будем ссылаться на него, как на динамический коэффициент диффузии (D-din). На фиг. 1 представлены экспериментально измеренные [27] значения действительной и мнимой частей параметра Фейбельмана ${{d}_{ \bot }}$ для золота в зависимости от длины волны. Можно видеть значительные изменения этого параметра в диапазоне от 400 до 600 нм.

Фиг. 1.

Зависимость действительной и мнимой частей параметра Фейбельмана ${{d}_{ \bot }}$ от длины волны.

При представлении численных результатов будут исследованы следующие модели.

• Классическая теория Максвелла – локальный случай (LRA).

• Теория GNOR с постоянным коэффициентом диффузии D.

• Модифицированная теория GNOR с динамическим коэффициентом диффузии (D-din).

• Модель с мезоскопическими граничными условиями, описываемыми параметрами Фейбельмана (SRF).

Будем исследовать задачу дифракции плоской волны р-поляризации на сфероидальной частице эквиобъемного диаметра D = 10 нм. Пусть частица располагается в воде с показателем преломления ${{n}_{e}} = 1.33$. Мы рассмотрим вытянутый сфероид с соотношением продольной и поперечной осей 1 : 3 (r = 3.0) и сплюснутый сфероид с соотношением осей 3 : 1 (r = 3.0). Во всех случаях плоская волна падает таким образом, чтобы вектор электрического поля был параллелен большей оси сфероида.

На фиг. 2 показано сечение экстинкции в зависимости от длины волны падающего излучения, рассчитанное в рамках классической теории Максвелла (локальный отклик) и модели SRF для двух различных значений параметра Фейбельмана. В данном случае результаты, полученные с помощью МДИ, сравнивались со значениями, которые даeт обобщенная неклассическая теория Ми [14] для сферической золотой частицы диаметром D = 10 нм. Отметим полное соответствие результатов МДИ и теории Ми. В то же время наблюдается существенное отличие графиков, соответствующих различным значениям параметра Фейбельмана: при ${{d}_{ \bot }}$ = +0.5 нм максимальное значение кривой SRF оказывается почти вдвое больше, чем у кривой LRA, при этом максимум смещается в сторону длинных волн, а при ${{d}_{ \bot }} = - 0.5$ нм, график SRF, наоборот, имеет сниженное значение максимума, который сдвигается в сторону коротких волн.

Фиг. 2.

Сечение экстинкции золотой сферы диаметром 10 нм, расположенной в ${\text{Si}}{{{\text{O}}}_{2}},{\text{\;}}$полученное в рамках моделей LRA и SRF методом дискретных источников и с помощью теории Ми.

На фиг. 3 представлены графики для сечения экстинкции (8) и коэффициента усиления (9), построенные для сферической частицы из золота диаметром 10 нм, находящейся в воде, и соответствующие различным моделям: LRA, GNOR, D-din и SRF. Ближе всего к локальному случаю (LRA) оказалась кривая SRF, у которой несколько уменьшен и сдвинут влево максимум. Кривые модели GNOR имеют максимальное значение, примерно на 40% меньше, чем у кривой LRA, причем учет динамического коэффициента диффузии корректирует максимальное значение в сторону увеличения по сравнению со стандартной моделью GNOR, использующей постоянное значение коэффициента. Наблюдается также небольшое смещение максимумов этих кривых в сторону коротких волн. Указанные особенности проявляются в большей степени на графиках, соответствующих коэффициенту усиления.

Фиг. 3.

Сечение экстинкции (a) и коэффициент усиления (б) золотой сферы диаметром 10 нм, расположенной в воде. Сравнение результатов, полученных в рамках моделей LRA, GNOR, D-din и SRF.

Фигура 4 демонстрирует сравнение аналогичных результатов, полученных в рамках моделей LRA, GNOR, D-din и SRF для сфероидальной частицы из золота, имеющей эквиобъемный диаметр 10 нм, соотношение поперечной и продольной осей 1 : 3 и расположенной в воде. Волна падает перпендикулярно продольной оси, т.е. угол падения ${{\theta }_{0}} = 90^\circ {\text{\;}}.$ Различия графиков, отмеченные на фиг. 3, проявляются в данном случае намного значительнее, причем в большей степени для коэффициента усиления. Максимальное значение кривой SRF оказывается почти в два раза меньше, чем в локальном случае (LRA), и сдвигается влево более чем на 25 нм. Снижение максимумов на графиках GNOR и D-din оказывается еще больше – более 4 раз по сравнению со значением LRA для коэффициента усиления, но сдвиг влево получается меньше – около 10 нм.

Фиг. 4.

Сечение экстинкции и коэффициент усиления для вытянутого золотого сфероида с соотношением осей 1: 3 и эквиобъемным диаметром 10 нм, расположенного в воде. Угол падения плоской волны – ${{{{\theta }}}_{0}} = 90^\circ .$ Сравнение результатов, полученных в рамках моделей LRA, GNOR, D-din и SRF.

На фиг. 5 показаны результаты сравнения моделей LRA, GNOR, D-din и SRF для сфероидальной частицы из золота с эквиобъемным диаметром 10 нм и соотношением поперечной и продольной осей 3 : 1. В данном случае волна падает перпендикулярно поперечной оси сфероида под углом ${{\theta }_{0}} = 0^\circ {\text{\;}}.$ В отличие от предыдущего случая, максимум кривой SRF снижается незначительно, но смещается теперь в другую сторону – на несколько нанометров в сторону длинных волн. Кривые GNOR и D-din демонстрируют значительное расхождение, которого не было ранее. Учет динамического коэффициента диффузии (D-din) приводит к существенному увеличению максимального значения по сравнению с традиционной моделью GNOR. При этом обе кривые, как и на фиг. 4, смещаются на несколько нанометров в сторону коротких волн относительно LRA.

Фиг. 5.

Сечение экстинкции и коэффициент усиления для сплюснутого золотого сфероида с соотношением осей 3 : 1 и эквиобъемным диаметром 10 нм, расположенного в воде. Угол падения плоской волны – ${{{{\theta }}}_{0}} = 0^\circ .$ Сравнение результатов, полученных в рамках моделей LRA, GNOR, D-din и SRF.

Рассмотрим теперь различия в результатах моделей LRA, GNOR, D-din и SRF на примере серебряной сфероидальной частицы с эквиобъемным диаметром 10 нм и соотношением поперечной и продольной осей 3 : 1 (фиг. 6). Частица размещена в воде, а плоская волна падает перпендикулярно поперечной оси сфероида под углом ${{\theta }_{0}} = 0^\circ {\text{\;}}.$ Следует отметить почти полное совпадение кривых SRF и LRA, что означает очень слабое влияние поверхностных эффектов в серебре. Напротив, кривые GNOR и D-din, демонстрирующие влияние объемных эффектов, ведут себя иначе. Максимум кривой GNOR оказывается в разы меньше, чем у остальных кривых, а учет динамического коэффициента диффузии (D-din) приводит к многократному росту максимальных значений. В то же время максимумы обеих кривых смещены влево относительно LRA, и больше всего сдвигается кривая GNOR – примерно на 12 нм.

Фиг. 6.

Сечение экстинкции и коэффициент усиления для сплюснутого серебряного сфероида с соотношением осей 3 : 1 и эквиобъемным диаметром 10 нм, расположенного в воде. Угол падения плоской волны – ${{{{\theta }}}_{0}} = 0^\circ $. Сравнение результатов, полученных в рамках моделей LRA, GNOR, D-din и SRF.

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

Сформулируем основные результаты статьи.

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

2. Сравнительный анализ объемных (GNOR) и поверхностных эффектов (SRF) позволил установить их существенные отличия. Если учет GNOR для благородных металлов всегда сопровождается снижением амплитуды и сдвигом в коротковолновую область, то проявление эффекта SRF существенно зависит от геометрии частицы, при этом он в значительной степени восстанавливает амплитуду ПР, сопровождая сдвигом, который может превышать 25 нм. Это различие особенно заметно при сравнительном анализе коэффициентов усиления поля на поверхности частицы.

3. Установлено, что различие между классическим случаем (LRA) и мезоскопическим (SRF) для серебра минимально. Это связано с тем, что для серебра параметр Фейбельмана отличен от нуля в узком диапазоне длин волн вблизи $\lambda = 250\,\,~{\text{нм}}$. В то же время для золота область изменения существенно сдвинута в видимый диапазон $\lambda = 450~\,\,{\text{нм}}$ [27].

Список литературы

  1. Barnes W.L., Dereux A., Ebbesen T.W. Surface plasmon subwavelength optics // Nature. 2003. V. 424. P. 824.

  2. Chon J.W.M., Iniewski K. Nanoplasmonics. Advanced Device Applications. CRC Press. 2018.

  3. Shi H., Zhu X., Zhang S. et al. Plasmonic metal nanostructures with extremely small features: new effects, fabrication and applications // Nanoscale Adv. 2021. V. 3. P. 4349.

  4. David C., Garcìa de Abajo F.J. Surface Plasmon Dependence on the Electron Density Profile at Metal Surfaces // ACS Nano. 2014. V. 8. № 9. P. 9558.

  5. Zhu W., Esteban R., Borisov A.G. et al. Quantum mechanical effects in plasmonic structures with subnanometre gaps // Nat. Commun. 2016. V. 7. 11495.

  6. Ullrich C.A. Time-Dependent Density-Functional Theory: Concepts and Applications. OUP Oxford. 2011.

  7. Sinha-Roy R., Garcìa-Gonźlez P., Weissker H.-C. et al. Classical and ab initio plasmonics meet at sub-nanometric noble metal rods // ACS Photonics. 2017. V. 4. № 6. P. 1484.

  8. Toscano G., Straubel J., Kwiatkowski A. et al. Resonance shifts and spill-out effects in self-consistent hydrodynamic nanoplasmonics // Nat. Commun. 2015. V. 6. № 1. P. 7132.

  9. Mortensen N.A., Raza S., Wubs M. et al. A generalized non-local optical response theory for plasmonic nanostructures // Nat. Commun. 2014. V. 5. 3809.

  10. Kupresak M., Zheng X., Gae V., Moshchalkov V.V. Appropriate nonlocal hydro- dynamic models for the characterization of deep-nanometer scale plasmonic scatterers // Adv. Theory Simul. 2019. V. 3. 1900172.

  11. Feibelman P.J. Surface electromagnetic fields // Prog. Surf. Sci. 1982. V. 12. 287.

  12. Deng H.-Y. A theory of electrodynamic response for bounded metals: Surface capacitive effects // Ann. Phys. 2020. V. 418. 168204.

  13. Yang Y., Zhu D., Yan W. et al. A general theoretical and experimental framework for nanoscale electromagnetism // Nature (London). 2019. V. 576. 248.

  14. Gonçalves P.A.D., Christensen T., Rivera N. et al. Plasmon–emitter interactions at the nanoscale // Nat. Commun. 2020. V. 11. 366.

  15. Stamatopoulou P.E., Tserkezis C. Finite-size and quantum effects in plasmonics: manifestations and theoretical modelling [Invited] // Optical Materials Express. 2022. V. 12. № 5. P. 1869.

  16. Mortensen N.A. Mesoscopic electrodynamics at metal surfaces (Review) // Nanophotonics 2021. V. 10. № 10. P. 2563.

  17. Yang F., Ding K. Transformation optics approach to mesoscopic plasmonics // Phys. Rev. B. 2022. V. 105. L121410.

  18. Mortensen N.A., Gonçalves P.A.D., Shuklin F.A. et al. Surface-response functions obtained from equilibrium electron-density profiles // Nanophotonics. 2021. V. 10. № 14. P. 3647.

  19. Еремин Ю.А., Свешников А.Г. Математическая модель учета эффекта нелокальности плазмонных структур на основе метода дискретных источников // Ж. вычисл. матем. и матем. физ. 2018. Т. 58. № 4. С. 586.

  20. Еремин Ю.А., Свешников А.Г. Квазиклассические модели квантовой наноплазмоники на основе метода Дискретных источников (обзор) // Журн. вычисл. матем.и матем. физ. 2021. Т. 61. № 4. С. 34.

  21. Doicu A., Eremin Yu., Wriedt T. Acoustic and Electromagnetic Scattering Analysis Using Discrete Sources. San Diego: Academic Press, 2000.

  22. Еремин Ю.А., Свешников А.Г. Математические модели задач нанооптики и биофотоники на основе метода дискретных источников // Ж. вычисл. матем. и матем. физ. 2007. Т. 47. № 2. С. 266.

  23. Бахвалов Н.С. Численные методы. М.: Наука, 1975.

  24. Колтон Д., Кресс Р. Методы интегральных уравнений в теории рассеяния. М.: Мир, 1987.

  25. Еремин Ю.А., Захаров Е.В. Аналитическое представление для интегрального поперечника рассеяния в рамках интегрофункционального метода Дискретных источников // Дифференц. ур-ния. 2022. Т. 58. № 8. С. 1073.

  26. Svendsen M.K., Wolff C., Jauho A.-P. et al. Role of diffusive surface scattering in nonlocal plasmonics // J. Phys.: Condens. Matter. 2020. V. 32. 395702.

  27. Echarri R.A., Gonçalves P.A.D., Tserkezis C. et al. Optical response of noble metal nanostructures: Quantum surface effects in crystallographic facets // Optica. 2021. V. 8. № 5. P. 710.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал