Доклады Российской академии наук. Математика, информатика, процессы управления, 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
Аннотация
На основе метода дискретных источников разработан и реализован строгий подход, позволяющий проводить анализ рассеяния светового излучения моделями резонаторов плазмонного нанолазера (SPASER). Данный подход позволяет учитывать все особенности граничной задачи для системы Максвелла, включая взаимодействие резонаторов с поверхностью призмы, а также эффект нелокального экранирования в рамках модели обобщенного нелокального отклика. Показано, что модель резонатора с диэлектрическим ядром и плазмонной оболочкой имеет существенные преимущества перед слоистой моделью с плазмонным ядром. Установлены условия, позволяющие получить усиление интенсивности поля на несколько порядков. Показано, что учет обобщенного нелокального отклика приводит к снижению интенсивности поля на 50%.
Одной из фундаментальных проблем квантовой наноплазмоники является разработка и реализация наноразмерных источников когерентного излучения [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. Тогда математическая постановка задачи рассеяния с учетом ЭН в слое может быть записана как
(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 .$(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} $Для построения полей внутри ${{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} $Здесь $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} $$k_{{nl}}^{2} = \frac{{{{\varepsilon }_{s}}(\omega )}}{{{{\eta }^{2}}}}.$
Тогда приближенное решение для Р-поляризации принимает вид
(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 как
Тогда (θ, φ)-компоненты диаграммы на единичной полусфере для 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} $Конкретные выражения для спектральных функций (8) приведены в [12].
Таким образом, определив амплитуды ДИ для рассеянного поля, можно легко вычислить компоненты диаграммы направленности (7) в верхнем полупространстве на единичной полусфере ${{\Omega }^{ + }} = \{ 0 \leqslant \theta \leqslant \pi {\text{/}}2;0 \leqslant \varphi \leqslant 2\pi \} $. Интенсивность рассеянного поля определяется как
Полное сечение рассеяния, представляющее собой суммарную интенсивность рассеянного поля в верхнее полупространство, будет
(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}}$, т.е.
На рис. 3 приведены результаты исследования коэффициента усиления для различных толщин золотой оболочки d = 5, 2 нм. Видно, что уменьшение толщины приводит к сдвигу ПР в длинноволновую область и существенному росту коэффициента усиления. Учет ЭН влечет снижение максимума почти на 50%.
Зададимся вопросом: что еще можно предложить для усиления поля. Рисунок 4 посвящен результатам анализа зависимости коэффициента усиления от угла падения плоской волны. Хорошо видно, что на длине волны ПР $\lambda = 635$ нм увеличение наклона падения волны приводит к дополнительному росту коэффициента усиления до 30%. При этом учет ЭН снова существенно снижает величину F.
В заключение отметим, что в данной работе разработан и реализован строгий метод, позволяющий проводить анализ усиления поля моделями резонаторов плазмонного нанолазера с учетом эффекта нелокального экранирования. Данный метод позволяет учитывать все особенности граничной задачи, включая взаимодействие слоистых резонаторов с поверхностью призмы, а также эффект нелокальности в рамках модели обобщенного нелокального отклика. Показано, что модель резонатора с диэлектрическим ядром и плазмонной оболочкой имеет существенные преимущества перед слоистой моделью с плазмонным ядром. Установлены условия, позволяющие получить усиление интенсивности поля на несколько порядков. Показано, что учет обобщенного нелокального отклика приводит к снижению интенсивности усиления поля на 50%.
Список литературы
Балыкин В.И. // Успехи физ. наук. 2018. Т. 188. № 9. С. 935–963.
Xu D., Xiong X., Wu L., et al. // Advances Opt. Photon. 2018. V. 10. № 4. P. 703–756.
Stockman M.I., Kneipp K., Bozhevolnyi S.I., et al. // J. Opt. 2018. V. 20. N043001.
Bergman D.J., Stockman M.I. // Phys. Rev. Lett. 2003. V. 90. N027402.
Noginov M.A, Zhu G., Belgrave A.M., et al. // Nature. 2009. V. 460. P. 1110–1113.
Premaratne M., Stockman M. // Advances Opt. Photon. 2017. V. 9. № 1. P. 79–128.
Mortensen N.A., Raza S., Wubs M., Søndergaard T., BozhevolnyiS.I. // Nature Commun. 2014. V. 5. P. 3809–3815.
Barbry M., Koval P., Marchesin F., et al. // Nano Lett. 2015. V. 15. N3410.
Wubs M., Mortensen A. // Quantum Plasmonics. S.I. Bozhevolnyi (eds.). Springer, Switzerland, 2017. P. 279–302.
Еремин Ю.А., Свешников А.Г. // ЖВМиМФ. 2018. Т. 58. № 4. С. 586–594.
Mortensen N.A., Raza S., Wubs M., et al. // Nat. Commun. 2014. № 5. N3809.
Еремин Ю.А., Свешников А.Г. // ДАН. 2017. Т. 477. № 2. P. 153–158.
Дмитриев В.И., Захаров Е.В. Метод интегральных уравнений в вычислительной электродинамике. М.: Макс Пресс, 2008.
Дополнительные материалы отсутствуют.
Инструменты
Доклады Российской академии наук. Математика, информатика, процессы управления