Доклады Российской академии наук. Математика, информатика, процессы управления, 2020, T. 490, № 1, стр. 24-28

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

Ю. А. Еремин 1*, А. Г. Свешников 1**

1 Московский государственный университет им. М.В. Ломоносова
Москва, Россия

* E-mail: eremin@cs.msu.ru
** E-mail: sveshnikov@phys.msu.ru

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

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

Аннотация

На основе метода дискретных источников разработан и реализован строгий подход, позволяющий проводить анализ рассеяния светового излучения моделями резонаторов плазмонного нанолазера (SPASER). Данный подход позволяет учитывать все особенности граничной задачи для системы Максвелла, включая взаимодействие резонаторов с поверхностью призмы, а также эффект нелокального экранирования в рамках модели обобщенного нелокального отклика. Показано, что модель резонатора с диэлектрическим ядром и плазмонной оболочкой имеет существенные преимущества перед слоистой моделью с плазмонным ядром. Установлены условия, позволяющие получить усиление интенсивности поля на несколько порядков. Показано, что учет обобщенного нелокального отклика приводит к снижению интенсивности поля на 50%.

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

Одной из фундаментальных проблем квантовой наноплазмоники является разработка и реализация наноразмерных источников когерентного излучения [1]. Идея состоит в том, чтобы использовать плазмонные поля вместо фотонных, используемых в обычных лазерах. Как известно, плазмонные поля позволяют преодолеть дифракционное ограничение на размер лазера. Плазмонный нанолазер получил название спасер (SPASER – Surface Plasmon Amplification by Simulated Emission of Radiation) [2, 3]. Его концепция впервые предложена Стокманом и Бергманом в 2003 г. [4]. В первой экспериментальной реализации спасера использовалась слоистая структура, состоящая из золотой наносферы со сферической оболочкой из прозрачного диэлектрика [5]. До настоящего времени ведутся активные исследования, посвященные разработке различных перспективных схем плазмонных нанолазеров [1, 3].

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

Поскольку размеры слоистых частиц составляют несколько десятков нанометров, то описание полей в рамках классической теории Максвелла оказывается недостаточным, так как начинают проявляться квантово-механические эффекты, такие как эффект нелокального экранирования [7]. Для описания подобных эффектов используются чисто квантовые подходы, например функциональная теория плотности во временной области [8], однако наиболее востребованными являются квазиклассические модели, описывающие квантовые эффекты, так как они позволяют правильно описывать поведение оптических характеристик частиц диаметром несколько десятков нанометров [9].

В настоящей работе метод дискретных источников (МДИ) [10] обобщается для исследования влияния эффекта нелокальности (ЭН) в рамках модели обобщенного нелокального отклика (GNOR) [11] на характеристики резонаторов плазмонного нанолазера. На основе компьютерной реализации установлено, что учет GNOR может снижать усиление поля резонатора на 50%.

Перейдем к математической постановке задачи. Пусть все пространство ${{\mathbb{R}}^{3}}$ разделено на два полупространства: активная среда – D0 (z > 0) и подложка –D1 (z < 0). Обозначим $\Sigma $ (z = 0) плоскую границу раздела. Пусть слоистая частица целиком располагается в верхнем полупространстве D0. Внутреннюю область будем обозначать ${{D}_{i}}$, область слоя – ${{D}_{s}}$, а внутреннюю и внешние поверхности – $\partial {{D}_{{i,s}}} \in {{C}^{{(2,\alpha )}}}$. Пусть $\{ {{{\mathbf{E}}}^{0}},{{{\mathbf{H}}}^{0}}\} $ – поле плоской электромагнитной волны, распространяющейся под углом $\pi - {{\theta }_{0}}$ к оси 0z. Тогда математическая постановка задачи рассеяния с учетом ЭН в слое может быть записана как

$\begin{gathered} {\text{rot}}{{{\mathbf{H}}}_{\zeta }} = jk{{\varepsilon }_{\zeta }}{{{\mathbf{E}}}_{\zeta }};\quad {\text{rot}}{{{\mathbf{E}}}_{\zeta }} = - jk{{\mu }_{\zeta }}{{{\mathbf{H}}}_{\zeta }}\quad {\text{в}}\quad {{D}_{\zeta }}, \\ \zeta = 0,1,i, \\ \end{gathered} $
$\begin{gathered} {\text{rot}}{{{\mathbf{H}}}_{s}} = jk({{\varepsilon }_{s}} + {{\eta }^{2}}\nabla {\text{div}}){{{\mathbf{E}}}_{s}}(M); \\ {\text{rot}}{{{\mathbf{E}}}_{s}} = - jk{{\mu }_{s}}{{{\mathbf{H}}}_{s}}\quad {\text{в}}\quad {\text{ }}{{D}_{s}}, \\ \end{gathered} $
${{{\mathbf{n}}}_{i}} \times \left( {{{{\mathbf{E}}}_{i}}\left( P \right) - {{{\mathbf{E}}}_{s}}\left( P \right)} \right) = 0,$
${{{\mathbf{n}}}_{i}} \times \left( {{{{\mathbf{H}}}_{i}}\left( P \right) - {{{\mathbf{H}}}_{s}}\left( P \right)} \right) = 0,\quad P \in {\text{ }}\partial {{D}_{i}}{\text{;}}$
${{{\mathbf{n}}}_{i}} \cdot \left( {{{\varepsilon }_{i}}{{{\mathbf{E}}}_{i}}(P) - {{\varepsilon }_{{nl}}}{{{\mathbf{E}}}_{s}}\left( P \right)} \right) = 0.$
${{{\mathbf{n}}}_{s}} \times \left( {{{{\mathbf{E}}}_{s}}\left( Q \right) - {{{\mathbf{E}}}_{0}}\left( Q \right)} \right) = {{{\mathbf{n}}}_{s}} \times {\mathbf{E}}_{0}^{0}\left( Q \right),$
$\begin{gathered} {{{\mathbf{n}}}_{s}} \times \left( {{{{\mathbf{H}}}_{s}}\left( Q \right) - {{{\mathbf{H}}}_{0}}\left( Q \right)} \right) = {{{\mathbf{n}}}_{s}} \times {\mathbf{H}}_{0}^{0}\left( Q \right), \\ Q \in {\text{ }}\partial {{D}_{s}}{\text{;}} \\ \end{gathered} $
${{{\mathbf{n}}}_{s}} \cdot \left( {{{\varepsilon }_{{nl}}}{{{\mathbf{E}}}_{s}}(Q) - {{\varepsilon }_{0}}{{{\mathbf{E}}}_{0}}\left( Q \right)} \right) = {{\varepsilon }_{0}}{{{\mathbf{n}}}_{s}} \cdot {\mathbf{E}}_{0}^{0}(Q).$
(1a)
$\begin{array}{*{20}{c}} {{{{\mathbf{e}}}_{z}} \times \left( {{{{\mathbf{E}}}_{0}}(\gamma ) - {{{\mathbf{E}}}_{1}}(\gamma )} \right) = 0,} \\ {{{{\mathbf{e}}}_{z}} \times \left( {{{{\mathbf{H}}}_{0}}(\gamma ) - {{{\mathbf{H}}}_{1}}(\gamma )} \right) = 0,} \end{array}\quad {\text{ }}\gamma \in \Sigma .$
$\begin{gathered} \mathop {\lim }\limits_{r \to \infty } {\text{ }}r \cdot \left( {{{{\mathbf{H}}}_{\zeta }} \times \frac{{\mathbf{r}}}{r} - {{{\mathbf{E}}}_{\zeta }}} \right) = 0,\quad r = \left| M \right| \to \infty , \\ \zeta = 0,1,\quad z \ne 0; \\ \end{gathered} $
(1б)
$\begin{gathered} \max \left( {\left| {{{{\mathbf{H}}}_{\zeta }}} \right|,\left| {{{{\mathbf{E}}}_{\zeta }}} \right|} \right) = O(\sqrt \rho ),\quad \rho = \sqrt {{{x}^{2}} + {{y}^{2}}} , \\ \rho \to \infty ,\quad z = \pm 0. \\ \end{gathered} $

Здесь $\{ {{{\mathbf{E}}}_{{0,1}}},{{{\mathbf{H}}}_{{0,1}}}\} $ – рассеянное поле в ${{D}_{{0,1}}}$, $\{ {{{\mathbf{E}}}_{{i,s}}},{{{\mathbf{H}}}_{{i,s}}}\} $ – поле в ${{D}_{{i,s}}}$, ${{{\mathbf{E}}}_{s}} = {\mathbf{E}}_{s}^{T} + {\mathbf{E}}_{s}^{L}$, ${\text{div}}{\mathbf{E}}_{s}^{T} = 0$, ${\text{rot}}{\mathbf{E}}_{s}^{L}$ = 0, ${{{\mathbf{n}}}_{{i,s}}}$ – единичные нормали к $\partial {\kern 1pt} {{D}_{{i,s}}}$, ${{{\mathbf{e}}}_{z}}$ – орт 0z, $\{ {\mathbf{E}}_{0}^{0},{\mathbf{H}}_{0}^{0}\} $ связано с внешним возбуждением и будет определено ниже, k = ω/c, а характеристики среды удовлетворяют условиям $\operatorname{Im} {{\varepsilon }_{{0,1,i}}},{{\mu }_{{0,1,i}}} = 0$, Imεs, ${{\mu }_{s}} \leqslant 0$, $\operatorname{Im} {{\varepsilon }_{{nl}}} = 0$. Временная зависимость была выбрана в виде $\exp \left\{ {j\omega {\text{ }}t} \right\}$. Параметры $\eta $ и ${{\varepsilon }_{{nl}}}$, относящиеся к продольному полю ${\mathbf{E}}_{s}^{L}$, определяются как ${{\varepsilon }_{{nl}}} = {{\varepsilon }_{s}} - \frac{{\omega _{p}^{2}}}{{j\gamma \omega - {{\omega }^{2}}}}$, η2 = $\frac{{{{\varepsilon }_{{nl}}}\left( {{{\beta }^{2}} + D\left( {\gamma + j\omega } \right)} \right)}}{{{{\omega }^{2}} - j\gamma \omega }}$. Здесь ${{\omega }_{p}}$ – плазменная частота металла, $\gamma $ – коэффициент затухания, β – гидродинамическая скорость в плазме связана со скоростью Ферми ${{{v}}_{F}}$ соотношением ${{\beta }^{2}} = \frac{3}{5}{v}_{F}^{2}$, D – коэффициент диффузии электронов [11].

Будем строить приближенное решение задачи (1), руководствуясь схемой метода дискретных источников [12]. Сначала решим задачу отражения и преломления поля плоской волны $\{ {{{\mathbf{E}}}^{0}},{{{\mathbf{H}}}^{0}}\} $ на $\Sigma $ и суммарное поле падающей и отраженной плоских волн в D0 обозначим как $\{ {\mathbf{E}}_{0}^{0},{\mathbf{H}}_{0}^{0}\} $. При построении приближенного решения задачи (1) для полей будем учитывать осевую симметрию и поляризацию внешнего возбуждения [11]. Суть такого подхода состоит в представлении для поля в каждой из областей ${{D}_{\zeta }}$, $\zeta = 0,1,i,s$, в виде линейной комбинации полей мультиполей, которые аналитически удовлетворяют модифицированной системе Максвелла, условиям излучения, а также условиям сопряжения для тангенциальных компонент полей на Σ. В основу представления полей в ${{D}_{{0,1}}}$ положим векторные потенциалы, элементами которых служат фурье-компоненты тензора Грина полупространства [13], записанные в цилиндрической системе координат как

(2)
$\begin{gathered} {\mathbf{A}}_{{mn}}^{{(e)e}} = \{ G_{m}^{e}(\xi ,z_{n}^{e})\cos (m + 1)\varphi ; \\ - G_{m}^{e}(\xi ,z_{n}^{e}) \times \sin (m + 1)\varphi ; \\ - g_{m}^{e}(\xi ,z_{n}^{e})\cos (m + 1)\varphi \} ,\quad e = 0,1; \\ \end{gathered} $
$\begin{gathered} {\mathbf{A}}_{{mn}}^{{(h)e}} = \{ G_{m}^{h}(\xi ,z_{n}^{e})\sin (m + 1)\varphi ; \\ G_{m}^{h}(\xi ,z_{n}^{e}\cos (m + 1)\varphi ; - g_{{m + 1}}^{h}(\xi ,z_{n}^{e})\sin (m + 1)\varphi \} , \\ {\mathbf{A}}_{{0n}}^{{(e)e}} = \{ 0;0;{\text{G}}_{0}^{h}(\xi ,z_{n}^{e})\} . \\ \end{gathered} $

Для построения полей внутри ${{D}_{{i,s}}}$ будут использоваться следующие потенциалы:

(3)
$\begin{gathered} {\mathbf{A}}_{{mn}}^{{(e)\nu }} = \{ Y_{m}^{\nu }(\xi ,z_{n}^{\nu })\cos (m + 1)\varphi ;\; - Y_{m}^{\nu }(\xi ,z_{n}^{\nu }) \times \\ \times \;\sin (m + 1)\varphi ;\;{\text{0}}\} ,\quad \nu = i,s \pm ; \\ \end{gathered} $
${\mathbf{A}}_{{mn}}^{{(h)\nu }} = \{ Y_{m}^{\nu }(\xi ,z_{n}^{\nu }){\text{sin}}(m + 1)\varphi ;\;Y_{m}^{\nu }(\xi ,z_{n}^{\nu }){\text{cos}}(m + 1)\varphi ;\;{\text{0}}\} ,$
${\mathbf{A}}_{n}^{{(e)\nu }} = \{ 0;\;0;\;{\text{ }}Y_{0}^{\nu }(\xi ,z_{n}^{\nu })\} $.

Здесь $Y_{m}^{i}(\xi ,z_{n}^{i}) = {{j}_{m}}({{k}_{i}}{{r}_{{\xi {\text{z}}_{n}^{i}}}}){{\left( {\frac{\rho }{{{{r}_{{\xi {\text{z}}_{n}^{i}}}}}}} \right)}^{m}}$, ${{j}_{m}}( \cdot )$ – сферическая функция Бесселя, $Y_{m}^{{s \pm }}(\xi ,z_{n}^{{s \pm }})$ = $h_{m}^{{(2,1)}}{\text{(}}{{k}_{{s \pm }}}{{r}_{{\xi {\text{z}}_{n}^{{s \pm }}}}}{\text{)}}{{\left( {\frac{\rho }{{{{r}_{{\xi {\text{z}}_{n}^{{s \pm }}}}}}}} \right)}^{m}}$, $h_{m}^{{(2,1)}}( \cdot )$ – сферические функция Ханкеля, $r_{{\xi {\text{z}}_{n}^{\nu }}}^{2}$ = = ${{\rho }^{2}} + {{(z - z_{n}^{\nu })}^{2}}$, $\xi = (\rho ,z){\text{ }}$, ki, s =  $k\sqrt {{{\varepsilon }_{{i,s}}}{{\mu }_{{i,s}}}} $, $z_{n}^{\nu }$ – координаты дискретных источников (ДИ). Для случая P-поляризации продольное поле строится на основе следующих скалярных потенциалов:

(4)
$\begin{gathered} \Psi _{{mn}}^{{s \pm }}(M) = h_{{m + 1}}^{{(2,1)}}({{k}_{{nl}}}{{R}_{{\xi z_{n}^{s}}}})\cos (m + 1)\varphi , \\ m = 0,1,\; \ldots ,\;M,\quad n = 1,2,\; \ldots ,\;{{N}_{{nl}}}; \\ \Psi _{n}^{{s \pm }}(M) = h_{0}^{{(2,1)}}({{k}_{{nl}}}{{R}_{{\xi z_{n}^{s}}}}), \\ \end{gathered} $
где волновое число определяется следующим образом [11]:

$k_{{nl}}^{2} = \frac{{{{\varepsilon }_{s}}(\omega )}}{{{{\eta }^{2}}}}.$

Тогда приближенное решение для Р-поляризации принимает вид

$\begin{gathered} {\mathbf{E}}_{\nu }^{{NT}} = \sum\limits_{m = 0}^M {{\text{ }}\sum\limits_{n = 1}^{N_{\nu }^{m}} {\left\{ {p_{{mn}}^{\nu }\frac{j}{{k{{\varepsilon }_{\nu }}{{\mu }_{\nu }}}}{\text{rot}}\,{\text{rot}}{\mathbf{A}}_{{mn}}^{{(e)\nu }}} \right.} + } \\ + \;\left. {q_{{mn}}^{\nu }\frac{1}{{{{\varepsilon }_{\nu }}}}{\text{rot}}{\mathbf{A}}_{{mn}}^{{(h)\nu }}} \right\} + \sum\limits_{n = 1}^{N_{{e,i}}^{0}} {r_{n}^{\nu }\frac{j}{{k{{\varepsilon }_{\nu }}{{\mu }_{\nu }}}}{\text{rot}}\,{\text{rot}}{\mathbf{A}}_{n}^{{(e)\nu }};} \\ \end{gathered} $
(5)
$\begin{gathered} {\mathbf{E}}_{{s \pm }}^{{NL}} = \sum\limits_{m = 0}^M {\sum\limits_{n = 1}^{N_{s}^{m}} {\bar {p}_{{mn}}^{{s \pm }}\nabla \Psi _{{mn}}^{{s \pm }}(M) + } } \sum\limits_{n = 1}^{N_{s}^{m}} {\bar {r}_{n}^{{s \pm }}\nabla \Psi _{n}^{{s \pm }}(M)} ; \\ {\mathbf{H}}_{\nu }^{N}{\text{ }} = \frac{j}{{k{{\mu }_{\nu }}}}{\text{rot}}{\mathbf{E}}_{\nu }^{N},\quad \nu = 0,1,i,s \pm . \\ \end{gathered} $

Совершенно аналогично строится решение для S-поляризации [12].

Заметим, что внутри плазмонной оболочки поле строится как сумма “уходящих” ($h_{{m + 1}}^{{(2)}}$) и “приходящих” ($h_{{m + 1}}^{{(1)}}$) волн, т.е.

(6)
$\begin{gathered} {\mathbf{E}}_{s}^{N} = {\mathbf{E}}_{{s + }}^{{NT}} + {\mathbf{E}}_{{s - }}^{{NT}} + {\mathbf{E}}_{{s + }}^{{NL}} + {\mathbf{E}}_{{s - }}^{{NL}},\quad {\text{div}}{\mathbf{E}}_{{s \pm }}^{{NT}} = 0, \\ {\text{rot}}{\mathbf{E}}_{{s \pm }}^{{NL}} = 0;\quad {\mathbf{H}}_{s}^{N}{\text{ }} = \frac{j}{{k{{\mu }_{s}}}}{\text{rot}}{\mathbf{E}}_{s}^{N} \\ \end{gathered} $.

Следует отметить, что построенные поля (5) и (6) удовлетворяют всем условиям граничной задачи (1), включая условия сопряжения на бесконечной поверхности подложки Σ, за исключением условий сопряжения на поверхностях ядра $\partial {\kern 1pt} {{D}_{i}}$ и оболочки $\partial {\kern 1pt} {{D}_{s}}$ резонатора. Именно эти условия используются для определения неизвестных амплитуд ДИ $\{ p_{{mn}}^{\nu },{\text{ }}q_{{mn}}^{\nu },{\text{ }}r_{n}^{\nu };{\text{ }}\bar {p}_{{mn}}^{{s \pm }},{\text{ }}r_{n}^{{s \pm }}\} $, $\nu = 0,1,i,s \pm $. Вычислительный алгоритм определения амплитуд ДИ построен аналогично [10].

Для вычисления характеристик рассеяния в дальней зоне используется диаграмма направленности рассеянного поля F, которая определяется в D0 как

$\begin{gathered} {{E}_{e}}(r){\text{/}}\left| {{{E}^{0}}(z = 0)} \right| = \frac{{\exp \{ - j{{k}_{0}}r\} }}{r}{\mathbf{F}}(\theta ,\varphi ) + O\left( {\frac{1}{{{{r}^{2}}}}} \right), \\ r \to \infty ,\quad z > 0. \\ \end{gathered} $

Тогда (θ, φ)-компоненты диаграммы на единичной полусфере для P-поляризации принимают вид

(7)
$\begin{gathered} F_{\theta }^{P}(\theta ,\varphi ) = j{{k}_{0}}\sum\limits_{m = 0}^M {\cos \left( {\left( {m + 1} \right)\varphi } \right){{{(j{{k}_{0}}\sin \theta )}}^{m}}} \times \\ \times \;\sum\limits_{n = 1}^{N_{0}^{m}} {\{ p_{{nm}}^{0}[\bar {G}_{n}^{e}\cos \theta + j{{k}_{0}}\bar {g}_{n}^{e}{{{\sin }}^{2}}\theta ]} + q_{{nm}}^{0}\bar {G}_{n}^{h}\} - \\ - \;{\text{ }}j\frac{{{{k}_{0}}}}{{{{\varepsilon }_{0}}}}\sin \theta \sum\limits_{n = 1}^{N_{0}^{0}} {r_{n}^{0}\bar {G}_{n}^{h}} , \\ \end{gathered} $
$\begin{gathered} F_{\varphi }^{P}(\theta ,\varphi ) = - j{{k}_{0}}\sum\limits_{m = 0}^M {\sin \left( {(m + 1)\varphi } \right){{{(j{{k}_{0}}\sin \theta )}}^{m}}} \times \\ \times \;\sum\limits_{n = 1}^{N_{0}^{m}} {\{ p_{{nm}}^{0}\bar {G}_{n}^{e}} + q_{{nm}}^{0}[\bar {G}_{n}^{h}\cos \theta + j{{k}_{0}}\bar {g}_{n}^{h}{{\sin }^{2}}\theta ]\} , \\ \end{gathered} $
где соответствующие спектральные функции $\bar {G}_{n}^{{e,h}}$, $\bar {g}_{n}^{h}$ представляются следующими соотношениями:

$\begin{gathered} \bar {G}_{n}^{{e,h}}(\theta ) = \exp \{ j{{k}_{0}}z_{n}^{e}\cos \theta \} + \\ + \;{\text{A}}_{{11}}^{{e,h}}({{k}_{0}}\sin \theta ) \cdot \exp \{ - j{{k}_{0}}z_{n}^{e}\cos \theta \} ,\quad z_{n}^{e} > 0; \\ \end{gathered} $
$\bar {g}_{n}^{{e,h}}(\theta ) = j{{k}_{0}}\cos \theta {v}_{{31}}^{{e,h}}({{k}_{0}}\sin \theta ,z = 0,z_{n}^{e}).$

Конкретные выражения для спектральных функций (8) приведены в [12].

Таким образом, определив амплитуды ДИ для рассеянного поля, можно легко вычислить компоненты диаграммы направленности (7) в верхнем полупространстве на единичной полусфере ${{\Omega }^{ + }} = \{ 0 \leqslant \theta \leqslant \pi {\text{/}}2;0 \leqslant \varphi \leqslant 2\pi \} $. Интенсивность рассеянного поля определяется как

$DS{{C}^{{P,S}}}({{\theta }_{0}},\theta ,\varphi ) = {\text{|}}F_{\theta }^{{P,S}}({{\theta }_{0}},\theta ,\varphi ){{{\text{|}}}^{2}} + \,{\text{|}}F_{\varphi }^{{P,S}}({{\theta }_{0}},\theta ,\varphi ){{{\text{|}}}^{2}}.$

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

(9)
$\sigma _{{}}^{{P,S}}({{\theta }_{0}}) = \int\limits_{{{\Omega }^{ + }}} {DSC_{{}}^{{P,S}}{\text{(}}{{\theta }_{0}},\theta ,\varphi {\text{)}}} d\omega .$

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

Перейдем к численным результатам. В качестве вещества частицы будем рассматривать SiO2 с индексом рефракции ${{n}_{i}} = 1.46$, а в качестве пленки – золото (Au), для которого соответствующие квантовые значения выбраны равными $\hbar {{\omega }_{p}} = 9.02$ эВ, $\hbar \gamma = 0.071$ эВ, ${{{v}}_{F}} = 1.39$ мкм/с, $D = 8.62 \times {{10}^{8}}$ мкм2/с [9]. Пусть частица располагается на стеклянной подложке ${{n}_{1}} = 1.52$ в активной среде R6G ${{n}_{0}} = 1.326$. Задавая длину волны внешнего возбуждения λ, вычисляя соответствующее значение $\omega $, легко определить значения нелокальных параметров ${{\varepsilon }_{{nl}}}$, ${{k}_{{nl}}}$. Пусть внутренний диаметр слоистой частицы D = 14 нм, а золотая оболочка имеет толщину d по аналогии с экспериментальной демонстрацией спасера [5].

На рис. 1 приведены результаты расчета сечения рассеяния (9) в зависимости от длины волны λ для слоистой сферической частицы и угла падения ${{\theta }_{0}} = 0^\circ $. Рассматриваются две различные слоистые частицы, одна с золотым ядром и диэлектрической оболочкой (${\text{Au@Si}}{{{\text{O}}}_{{\text{2}}}}$) и другая – с диэлектрическим ядром и золотой оболочкой (${\text{Si}}{{{\text{O}}}_{{\text{2}}}}{\text{@Au}}$). Толщина оболочки в обоих случаях $d = 10$ нм. Из рис. 1 видно, что учет ЭН приводит к уменьшению амплитуды плазмонного резонанса (ПР) и к небольшому сдвигу его в коротковолновую область. Следует отметить, что модель ${\text{Si}}{{{\text{O}}}_{{\text{2}}}}{\text{@Au}}$ оказывается на два порядка более эффективной, чем использованная в экспериментах [5] модель ${\text{Au@Si}}{{{\text{O}}}_{{\text{2}}}}$. Аналогичные результаты, относящиеся к коэффициенту усиления поля на внешней границе оболочки $\partial {\kern 1pt} {{D}_{s}}$, т.е.

$F = \frac{{\int\limits_{\partial {{D}_{s}}} {{\text{|}}{\mathbf{E}}_{0}^{N} + {\mathbf{E}}_{0}^{0}{{{\text{|}}}^{2}}} d\sigma }}{{\int\limits_{\partial {{D}_{s}}} {{\text{|}}{\mathbf{E}}_{0}^{0}{\text{|}}} d\sigma }},$
изображены на рис. 2. Как явствует из рис. 2, модель ${\text{Si}}{{{\text{O}}}_{{\text{2}}}}{\text{@Au}}$ снова оказывается в десятки раз эффективней. Следует отметить, что для модели ${\text{Au@Si}}{{{\text{O}}}_{{\text{2}}}}$ с учетом ЭН усиление составляло всего 1%, т.е. F = 1.01. В дальнейшем будем при моделировании использовать исключительно модель ${\text{Si}}{{{\text{O}}}_{{\text{2}}}}{\text{@Au}}$.

Рис. 1.

SCS: (9) в зависимости от длины волны λ для слоистой частицы D = 14 нм, d = 10 нм при угле падения плоской волны ${{\theta }_{0}} = {{0}^{ \circ }}$. Модель ${\text{Au@Si}}{{{\text{O}}}_{{\text{2}}}}$: кривая 1 – без учета ЭН (LRA); кривая 2 – с учетом ЭН: GNOR, модель ${\text{Si}}{{{\text{O}}}_{{\text{2}}}}{\text{@Au}}$: 3 – LRA, 4 – GNOR.

Рис. 2.

Коэффициент усиления F в зависимости от длины волны λ для слоистой частицы. Модель ${\text{Au@Si}}{{{\text{O}}}_{{\text{2}}}}$: кривая 1 – LRA, кривая 2 – GNOR, модель ${\text{Si}}{{{\text{O}}}_{{\text{2}}}}{\text{@Au}}$: 3 – LRA, 4 – GNOR.

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

Рис. 3.

F-модель ${\text{Si}}{{{\text{O}}}_{{\text{2}}}}{\text{@Au}}$: кривая 1d = 5 нм, LRA; 2 – GNOR; 3d = 2 нм, LRA; 4 – GNOR.

Зададимся вопросом: что еще можно предложить для усиления поля. Рисунок 4 посвящен результатам анализа зависимости коэффициента усиления от угла падения плоской волны. Хорошо видно, что на длине волны ПР $\lambda = 635$ нм увеличение наклона падения волны приводит к дополнительному росту коэффициента усиления до 30%. При этом учет ЭН снова существенно снижает величину F.

Рис. 4.

F в зависимости от угла падения ${{\theta }_{0}}$ при $\lambda = 635$ нм и d = 2 нм, модель ${\text{Si}}{{{\text{O}}}_{{\text{2}}}}{\text{@Au}}$, кривая 1 – LRA; 2 – GNOR.

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

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

  1. Балыкин В.И. // Успехи физ. наук. 2018. Т. 188. № 9. С. 935–963.

  2. Xu D., Xiong X., Wu L., et al. // Advances Opt. Photon. 2018. V. 10. № 4. P. 703–756.

  3. Stockman M.I., Kneipp K., Bozhevolnyi S.I., et al. // J. Opt. 2018. V. 20. N043001.

  4. Bergman D.J., Stockman M.I. // Phys. Rev. Lett. 2003. V. 90. N027402.

  5. Noginov M.A, Zhu G., Belgrave A.M., et al. // Nature. 2009. V. 460. P. 1110–1113.

  6. Premaratne M., Stockman M. // Advances Opt. Photon. 2017. V. 9. № 1. P. 79–128.

  7. Mortensen N.A., Raza S., Wubs M., Søndergaard T., BozhevolnyiS.I. // Nature Commun. 2014. V. 5. P. 3809–3815.

  8. Barbry M., Koval P., Marchesin F., et al. // Nano Lett. 2015. V. 15. N3410.

  9. Wubs M., Mortensen A. // Quantum Plasmonics. S.I. Bozhevolnyi (eds.). Springer, Switzerland, 2017. P. 279–302.

  10. Еремин Ю.А., Свешников А.Г. // ЖВМиМФ. 2018. Т. 58. № 4. С. 586–594.

  11. Mortensen N.A., Raza S., Wubs M., et al. // Nat. Commun. 2014. № 5. N3809.

  12. Еремин Ю.А., Свешников А.Г. // ДАН. 2017. Т. 477. № 2. P. 153–158.

  13. Дмитриев В.И., Захаров Е.В. Метод интегральных уравнений в вычислительной электродинамике. М.: Макс Пресс, 2008.

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

Инструменты

Доклады Российской академии наук. Математика, информатика, процессы управления