CC BY
20
УДК 517.958:535.4
А.В. Барышев1, Ю.А. Еремин2
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ОПТИЧЕСКОЙ АНТЕННЫ НА ОСНОВЕ МЕТОДА ДИСКРЕТНЫХ ИСТОЧНИКОВ*
Рассматривается оптическая антенна, представляющая собой линейный кластер нано-размерных частиц, заключенных внутри металлической пленки на подложке. На основе метода дискретных источников разработана и реализована математическая модель оптической антенны. Численные результаты демонстрируют возможность управления как шириной лепестка, так и направленностью излучения антенны.
Ключевые слова: метод дискретных источников, оптическая антенна, наночастицы, плазмонный резонанс.
1. Введение. Оптические наноразмерные антенны имеют широкий спектр применений в лазерной оптике, спектроскопии, биосенсорике и оптических микроскопах сверхразрешения [1, 2]. Как правило, оптические антенны представляют собой конструкцию, состоящую из совокупности нано-размерных частиц, локализованных на поверхности подложки или расположенных в тонкой пленке, нанесенной на поверхность подложки [2]. Интерес к подобным устройствам обусловлен тем, что они усиливают взаимодействие света с веществом, позволяя таким образом повысить эффективность таких оптоэлектронных устройств, как светодиоды, локальные биосенсоры и переключатели излучения. В настоящее время в процессе математического моделирования функциональных свойств оптических антенн широко применяется метод конечных разностей во временной области [3]. Он позволяет проводить моделирование сложных систем, однако не учитывает в полной мере взаимодействие между частицами и бесконечными границами раздела сред подложка-пленка и пленка-воздух, что приводит к существенным погрешностям при расчетах. В настоящей работе на основе метода дискретных источников (МДИ) [4] разработана и реализована математическая модель оптической антенны, представляющей собой совокупность наноразмерных частиц, расположенных в тонкой металлической пленке, нанесенной на поверхность прозрачной подложки. Приведены численные результаты, иллюстрирующие возможность управления прошедшим световым излучением.
2. Математическая постановка задачи рассеяния. Пусть все пространство разделено параллельными плоскостями Н0 (г = 0) и Н1 (г = й) на три части: воздух — область 1)0 (г > с1), проводящая пленка — Df (й > г > 0) и стеклянная призма — 1?1 (г < 0). В качестве внешнего возбуждения будем рассматривать поле линейно поляризованной плоской волны {Е°,Н0}, падающей из призмы под углом 01 относительно нормали на плоскую границу Н1 раздела призма-пленка. Пленка занимает пространство между плоскостями Н0 и Н1, при этом плоскость Н0 отделяет пленку от воздуха; N частиц с гладкой границей произвольной формы, внутреннюю область которых обозначим /), (. гп = г 1,... расположены целиком внутри пленки. Полагаем, что области /), ( являются односвязными. Введем прямоугольную систему координат с началом на плоскости Н1 и осью Ог, перпендикулярной плоскости Н1 и направленной в сторону области 1)0. Обозначим границы частиц через ¿Шгп, гп = ¿1,... ,гдг- Тогда математическая постановка задачи будет иметь вид
ml Н^ jk:\Ec. TotE^ =—jkfi^H^ в D^, С = 0, /, 1, ..., ipf,
щп X (Ein(p)-E/(p)) = 0, щп X (Hin(p) -H/(p)) = 0, pedDi.
x (Eo(p) — E/(p)) = 0, ez x (H0(p) — H/(p)) = 0, p € H0,
егх(Е/(р)-Е1(р)) = 0, егх(Н/(р)-Н1(р)) = 0, p G Hi,
1 Физический факультет МГУ, acn., e-mail: eclecticaQrambler.ru
2 Факультет ВМК МГУ, вед. науч. сотр., д.ф.-м.н., e-mail: ereminQcs.msu.su
* Работа выполнена при финансовой поддержке РФФИ, проект № 09-01-00318.
(|Е/|,|Н/|) = о (ехр {— |Im р}), р = \J х2 + у2 ^ оо.
Здесь — полные, а — рассеянные поля в соответствующих областях Д^, ez —
орт z декартовой системы координат, — нормали к поверхностям dDin € С'^2,а\ in = ii,..., ¿jv-Мы также полагаем, что параметры сред удовлетворяют условиям 1т е^, ^ 0, Im£f,Pf < 0 (временная зависимость выбрана в виде exp{jwi}). В этом случае граничная задача (1) имеет единственное решение [5]. Отметим, что в областях Dij полное поле включает в себя падающую и отраженную плоские волны, а в области D0 — преломленную волну.
Прежде чем строить приближенное решение для рассеянных полей, решим задачу дифракции поля плоской волны {Е°,Н0} на слоистой структуре воздух-пленка-подложка. Это решение можно записать в явном виде [6]. Полученное поле обозначим через {Е°, £ = 1, /, 0. Теперь будем строить приближенное решение задачи (1) для рассеянного поля {Е^, Н^} в Д^, £ = 1, /, 0, и полного поля в Dj на основе метода дискретных источников (МДИ) [4].
3. Метод дискретных источников. Суть МДИ заключается в представлении поля в виде конечной линейной комбинации полей дискретных источников, которая удовлетворяет системе уравнений Максвелла в областях />i. /-.п./,.....v. условиям на бесконечности для рассеянного в -Di,/,o поля,
а также условиям сопряжения для тангенциальных компонент полей на Н0д. В результате чего решение граничной задачи рассеяния (1) сводится к аппроксимации поля внешнего возбуждения на поверхностях частиц dDin, гп = гi,... ,ijv, полями дискретных источников. Неизвестные амплитуды ДИ определяются только из условий сопряжения на поверхностях частиц dDin, которые принимают вид
nin х (Ein - Щ) = nin х Е°, nin х (Hjn - Щ) = nin х Н°, (2)
где ¡„ ¡\.....i}v, in = iъ • • •; Ьv-
Приближенное решение будем строить для рассеивателей произвольной формы. В основу для представления рассеянного поля положим систему электрических диполей, удовлетворяющих условиям сопряжения для тангенциальных компонент полей на плоскостях Нод, что позволяет учесть аналитически всевозможные переотражения полей между частицами и плоскостями раздела сред [4]. Для g-й частицы, q = 1,..., N, диполи локализованы в точках Mq и располагаются всюду плотно на некоторой вспомогательной поверхности Sq € С2'", {Мп}п=1 = Sq. Поверхность Sq находится строго внутри Diq. В данном случае структура полей будет определяться тензором Грина слоистой среды [7]:
_ /Gii 0 0 \
G(M,M0)= 0 Си 0 .
\дд/дхм дд/дум G33j
При этом компоненты тензора могут быть записаны в виде интегральных представлений Зоммер-фельда
ос
Gu(M,MQ) = J JQ(\r)v11(\,z,zQ)\d\, о
ос
G33(M,MQ) = J JQ(\r)v33(\,z,zQ)\d\, о
оо
g(M,M0) = J JQ(\r)v31(\,z,zQ)\d\.
о
Здесь R2M M = r2 + (z — zq)2, r2 = p2 + pfj — 2ppQcos((p — (po), Jq( ) — цилиндрическая функция Бесселя, (p0, (p0, zq) — цилиндрические координаты точки M0. Для спектральных функций г/ц(z, zn, A), v33{z,t zfi, А), г/31 (z, zn, А), обеспечивающих выполнение условий сопряжения на границах Н0д слоистой
среды, справедливы следующие представления:
(он ехр{^т?о(.г — (¿)), г ^
ехр{-г}} - гп\}/% + /Зц ехр{-г)}{й - г)} + хп ехрй ^ г ^ О, <5п ехр{?712:}, г < О,
{«33 ехр{-щ(х-й)}, г ^ ё,
ехр{-г}} - гп\}/% + /З3з ехр{-г)}{й - г)} + хзз ехр-^/.г}, й ^ г ^ О, 5ззехр{г?1г}, г < О,
(«31 ехр{^т?0(.г - (¿)), х ^ й,
/З31 ехр{^т?/((^ - г)} + хз1 ехр{-77/2:}, й ^ г ^ О,
<531 ехр{г?1г}, г < 0.
Здесь г]^ = А2 — Щ, Щ = к2еКоэффициенты «ц = ац(А, /Зц = /Зц(А, (¿), хи =
= Хи(А, (¿) определяются из условий сопряжения, которым удовлетворяют спектральные функции, на границах г = {0, й} [7]:
ы = 1
1 ди
11
-^31
/х дг
= 0,
= 0, [г/33] =
1 ду31'
е/х
1 ди-
зз
е
1
= 0,
ер,
г^и,
где [•] — скачок переменной при переходе через поверхность.
Для д-й частицы, д = 1,..., Ж, полагаем, что в каждой точке М® располагаются три линейно независимых диполя, ориентированных в соответствии с цилиндрической системой координат. Тогда векторные потенциалы д-й частицы, связанные с представлением рассеянного поля и соответствующие каждому такому источнику, можно представить как
А^((г,е) = Сц(М, М%)соз((р - (рп)ер - Сц(М, М«) зт((р - (рп)е^ +
К ' п> соз(у> - срп)---^ ' п> 81п(у - срп)
ар р оср
А2(9'е) = Сц(М,М«)8тО - <рп)ер + Сц(М, М«) сов(<р - <рп)<
' " зт ^ - срп) + - -соз -стр Р О^Р
= С33(М,М«)ег.
Векторные потенциалы д-й частицы, д = 1,..., Ж, связанные с представлением полного поля внутри частицы, могут быть записаны следующим образом:
Ап9,г) = ММ, М«) С0в((р - срп)ер - ]о{М, М«) 8 т(ср - срп^ КШ) = ММ, М«) зт((р - (рп)ер + ММ, М«) со8(у> -
Здесь = jo(kiRмм%,)■> гДе ^о(') — сферическая функция Бесселя нулевого порядка.
Тогда для приближенного решения граничной задачи (1) справедливо следующее представление:
n ме
Е
(М) = А^'е)(М),
д=1п=1¿=1 А^ 3
Н 1'Н{М) = гot Е(М), М б Д
3
кр
е
(3)
/С/Хг
п=1 г=1
м € /),„.
где Е;"' д . Н;"'д — рассеянное системой частиц поле во внешней области /),. а Ед . — полное поле внутри д-й частицы, Ые, Щ — количество дискретных источников, выбираемых на вспомогательных поверхностях Бе, ш представляющих собой поле в соответствующих областях.
Построенное таким образом приближенное решение удовлетворяет всем условиям граничной задачи (1), за исключением условий сопряжения на поверхностях частиц. Удовлетворяя условиям (2) в некоторой норме, определим р^. Заметим, что с помощью единого набора амплитуд мы строим представление для рассеянного поля сразу во всех областях -С>1,/,о. Имеет место следующее утверждение.
Теорема. Пусть диполи локализованы в точках М®, д = 1,..., Ж, и располагаются всюду плотно на вспомогательной нерезонансной для данного волнового числа кй поверхности Б9 € С^2,а\
{Мп}п=1 = в9, при этом поверхность Б9 располагается строго внутри . Тогда приближенное решение вида (3) сходится равномерно к точному решению задачи (1) в любой конечной области полупространств £>од при Ме, Щ ^ оо.
4. Численный алгоритм. Как было отмечено выше, приближенное решение удовлетворяет большинству условий задачи (1), а определение амплитуд дискретных источников производится путем удовлетворения граничного условия (2) в норме Ь2(дО), где дБ = ("^¿Ш^. Для этого нам сначала потребуется построить поле внешнего возбуждения внутри пленки для Р/5-поляризованного
возбуждения плоской волны. Это поле может быть записано как
Е° = Е+ + и^Щ, Ш} = + wf's^ij,
где для Р-поляризации
Е^ = (ц^cosвfex + smвfez)'y±, Н^ = ^п/е|/7±
и для 5-поляризации
= соя вfex + smвfez)'y±, Е^ = е|/7±, ^ = exp{—jkf(xs'mвf ± гсояв^}, Щ = п^ = £ = 1,/, О,
и;1 = /(1 — е )-, и)г' = е ,
р с р с
^а'/З ' — соответствующие коэффициенты отражения и преломления на плоской границе, разделяющей /.),, и / ) >• [6].
Для вычисления неизвестных амплитуд дискретных источников применяются различные схемы. Установлено, что наиболее стабильные результаты получаются при псевдообращении переопределенной системы линейных уравнений, получаемой путем применения обобщенного метода коллокаций, суть которого заключается в том, что на поверхности каждой частицы мы выбираем набор точек коллокаций {Р^=1 € дО^, равномерно покрывающих эту поверхность [4]. Выполняя эту операцию для каждой поверхности частицы, а также задавая на вспомогательных поверхностях внутри частицы наборы дискретных источников и представляющих рассеянное поле вне частицы
и полное поле внутри частицы, мы получим переопределенную систему линейных уравнений размерности 11. N х ЗЖ(Же + N1). Затем неизвестные амплитуды мы определяем путем псевдообращения. Подобная схема позволяет одновременно вычислять характеристики рассеяния для обеих поляризаций и всего набора углов падения плоской волны.
Во всех наших численных экспериментах мы рассматривали линейный кластер сферических частиц, обладающий симметрией относительно плоскости '/.У. При этом дискретные источники располагались на вспомогательных сферах, радиус которых был вдвое меньше радиуса рассматриваемых сферических частиц, а численные значения координат и количество источников, представляющих собой рассеянное поле вне частиц и полное поле внутри частиц, совпадали.
Для вычисления характеристик рассеянного поля на бесконечности необходимо вычислить диаграмму рассеяния. В области 1)0 она определяется следующим образом:
Е%(М)/|Е°(г = 0)| = еХр(^оГ)Р02,у) + о(г~1), О г = \М\ ос.
Для получения конкретного вида диаграммы рассеяния для приближенного решения (3) достаточно использовать асимптотические представления для интегралов Зоммерфельда. Тогда в, ¡^-компоненты
диаграммы рассеяния записываются следующим образом:
N ne
F0(в, <р) = jko cosве^™в Y. ЕiPntfS + Pne2ÍÍ + Р%т,
q= 1 п=1
N ne
Fv(0, <р) = Jko COS Rodeóse J- ¿{рв-е/Р + p«'^},
g=l n=l
где
/g = {«иsiní?, zn, d) cos в + jkoahlko siní?, zn, d)sm20} cos(cp — (pn), fg = {«li siní?, zn,d) cos в + jkoah(ko siní?, zn, d)sm20} sin(9? — (pn), /I = -afjfco sin zm d) sin0, J? = ^afi(&osm0, zn,d) sin(9P - (pn), f$ = «и(^o sin0,zn,d) cos(cp - (pn).
Стоит также отметить, что зеркальная симметрия линейного кластера позволяет существенно сократить время вычислений. Это достигается за счет того, что матрица в левой части переопределенной системы линейных уравнений, из которой мы определяем неизвестные амплитуды дискретных источников, приводится к клеточно-циркулянтному виду. Применение клеточно-циркулянтного алгоритма позволило нам сократить время вычисления более чем в два раза в основном за счет того, что нам пришлось считать в два раза меньше интегралов Вейля-Зоммерфельда, на вычисление которых тратится основное время расчета.
5. Обсуждение результатов. Нас интересует интенсивность рассеянного поля в области D0, которая имеет вид
I(e^) = \F0(e^)\2 + \Fv(e^)\2.
В качестве внешнего возбуждения во всех численных экспериментах будем рассматривать Р-поляри-зованную плоскую волну единичной амплитуды с длиной волны А = 532 нм, волновой вектор которой лежит в плоскости (р = 0. Волна распространяется из призмы, сделанной из стекла с высоким индексом рефракции п\ = 1.904, Р-поляризация плоской волны была выбрана, поскольку в этом случае наблюдается эффект экстремального просачивания энергии [8]. Пусть во — преломленный угол, под которым волна проходит в D0. В этом случае соотношение Снеллиуса имеет вид smdo = (пг/щ) sin#i. Поскольку показатель преломления стекла больше показателя преломления воздуха, то при изменении угла падения 9i в диапазоне от 0 до (тг/2) наступает момент, начиная с которого |sin0o| > 1. Тогда
cosí?o = —j\/sin2í?o — 1 и 7+ = exp {^jfcoa;sin0o} exp — lj. Таким образом, в верхнем
полупространстве возникает неизлучающая волна, которая распространяется вдоль границы раздела (оси ОХ) и экспоненциально затухает по нормали к границе (вдоль оси OZ). При длине волны А = 532 нм для стекла с высоким индексом рефракции п\ = 1.904 критический угол вс = arcsin(no/ni), за которым располагается область неизлучающих волн, оказывается равным вс и 31.68.
Во всех численных экспериментах в качестве рассеивателей будем рассматривать кремниевые (Si) сферы диметром D = 40 нм, показатель преломления которых щ = (4.15; ^0.047) и которые целиком расположены внутри серебряной пленки (Ag) толщины d = 45 нм, имеющей показатель преломления щ = (0.14;—3.19). Толщина пленки выбрана из условия возникновения плазмонного резонанса [8]. Рассеяную интенсивность мы будем анализировать в плоскости падения плоской волны (р = 0: 1(6) = 1(6,ц> = 0).
На рис. 1 приведена интенсивность рассеяния для одной, двух и четырех сфер при нормальном падении волны, а также интенсивность для четырех сфер в случае, когда волна падает под углом 01 = 30°. Расстояние между соседними сферами равно 20 нм (расстояние между их центрами 1с = 60 нм). Видно, что с увеличением числа сфер интенсивность существенно возрастает.
На рис. 2 показано сравнение интенсивности для четырех сфер в случае нормального падения в зависимости от расстояния между частицами (S = 20,40,80 нм). Также приведена интенсивность для четырех сфер в случае наклонного падения под углом в = 30° при расстоянии между частицами, равном 80 нм. Можно заметить, что максимальная амплитуда при нормальном падении наблюдается
2
/(0), мкм
Рис. 1. Интенсивность рассеяния 1(0) в зависимости от количества сфер: 1 — одна сфера, 0\ =0°; 2 — две сферы, 8 = 20 нм, 0\ =0°; 3 — четыре сферы, 8 = 20 нм, 9\ = 0°; 4 — четыре сферы, 8 = 20 нм, 01 = 30°
2
/(0), МКМ
Рис. 2. Интенсивность рассеяния на четырех сферах 1(0) в зависимости от расстояния между сферами: 1 — 8 = 20 нм, 0 = 0°; 2 — 8 = 40 нм, 0 = 0°; 3 — 6 = 80 нм, 0 = 0°; 4 — 6 = 80 нм, 0 = 30°
при наименьшем расстоянии между частицами в 20 нм. В то же время интенсивность имеет более узкую направленность в случае, когда расстояние между частицами равно 80 нм. В случае наклонного падения при расстоянии между частицами в 80 нм интенсивность имеет меньшую амплитуду, чем в случае нормального падения, при этом возникает также небольшой побочный максимум.
Как известно [9], для того чтобы получить диаграмму направленности линейной решетки, необязательно считать взаимодействие между всеми элементами этой решетки. Можно посчитать диаграмму направленности одного или нескольких элементов решетки, а затем численные значения полученной диаграммы домножить на соответствующий множитель решетки. На рис. 3 изображено сравнение ин-
2
/(0), МКМ
Рис. 3. Сравнение интенсивности рассеяния на четырех сферах 1(0) ч полученной с учетом всех взаимодействий (кривые 1 — 6 = 40 нм, 3 — 8 = 80 нм) и умножением интенсивности рассеяния на двух сферах на множитель решетки (кривые 2 — 8 = 40 нм, 4 — 8 = 80 нм)
тенсивности для четырех сфер при нормальном падении в зависимости от расстояния между сферами (5 = 40, 80 нм). При этом в двух случаях численные значения интенсивности получаются путем учета всевозможных взаимодействий между четырьмя сферами, а в двух других — путем расчета диаграммы направленности для двух сфер с последующим умножением на множитель решетки из четырех элементов. Как мы можем видеть, при расстоянии между частицами в 80 нм наблюдается хорошее соответствие между интенсивностью, рассчитанной с учетом всех взаимодействий и с применением множителя решетки. Воспользуемся полученным результатом и рассчитаем диаграмму направленности для оптической антенны, представляющей собой линейный кластер из двадцати сфер, где расстояние между соседними сферами равно 80 нм. При этом в качестве базы мы будем использовать диаграмму направленности для четырех сфер с учетом всех взаимодействий, после чего умножим численные значения полученной диаграммы на множитель решетки из двадцати элементов.
2
1(0), мкм
90
Рис. 4. Интенсивность рассеяния на линейном кластере из двадцати сфер 1(0) в зависимости от угла падения волны 01 при расстоянии между сферами, равном 80 нм: 1 — 0 = 0°; 2 — 0 = 10°; 3 — 0 = 20°; 4 — 0 = 33.7° (плазмонный резонанс)
На рис. 4 приведено сравнение интенсивности для линейного кластера из двадцати сфер в зависимости от угла падения волны (в± = 0, 10, 20, 33.7°). Следует отметить, что по сравнению со случаем четырех сфер интенсивность для двадцати сфер становится более узконаправленной. При угле падения в\ = 33.7° наблюдается экстремальное просачивание энергии [8].
Как явствует из проведенного анализа, возможно управлять как шириной луча, так и его направлением, варьируя число сфер и угол падающей плоской волны в±.
СПИСОК ЛИТЕРАТУРЫ
1. В ей к о В. П., Вознесенский Н. В., Воронин Ю.М. и др. Лазерная технология формирования оптических антенн для ближнепольных микроскопов и исследование их характеристик // Изв. РАН. Сер. Физическая. 63. № 11. С. 1954-1963.
2. Bharadwaj Р., Deutsch В., Novotny L. Optical antennas // Advanced in Optics and Photonics. 2009. 1. P. 438-483.
3. Capoglu I. R., Smith G.S. A direct time-domain FDTD near-field-to-far-field transform in the presence of an infinite grounded dielectic slab // IEEE Trans. Antennas Propag. 2006. 52. N 12. P. 3805-3814.
4. Еремин Ю.А., Свешников А.Г. Математические модели задач нанооптики и биофотоники на основе метода дискретных источников // ЖВМиМФ. 2007. 47. № 2. С. 266-284.
5. Захаров Е. В. О единственности и существовании решений интегральных уравнений электродинамики неоднородных сред // Вычислительные методы и программирование. Сер. 24. М.: Изд-во МГУ, 1975. С. 37-42.
6. Chew W. С. Waves and Fields in Inhomogeneoues Media. N.Y.: IEEE Press, 1995.
7. Дмитриев В. И. Поля в слоистых средах. М.: Изд-во МГУ, 1963.
8. Гришина Н.В., Еремин Ю.А., Свешников А. Г. Математическая модель слоистой структуры с на-норазмерным отверстием // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 2008. № 4. С. 11-16.
9. Фельд Я.Н., Бененсон Л. С. Основы теории антенн. М.: Дрофа, 2007.
Поступила в редакцию 17.05.10
MATHEMATICAL MODEL OF AN OPTICAL ANTENNA BASED ON THE DISCRETE SOURCES METHOD
Baryshev A. V., Eremin Yu. A.
Optical antenna consisting of a linear cluster of nanodimensional particles located inside a metal film deposited on a substrate is considered. Based on Discrete Sources Method mathematical model of the optical antenna has been elaborated and realized. Numerical results demonstrate an ability to handle both width and direction of the antenna's radiation.
Keywords: Discrete Sources Method, optical antenna, nanoparticles, plasmonic resonance.
