Метод простых итераций с коррекцией сходимости и метод минимальных невязок в задачах плазмоники

Обложка

Цитировать

Полный текст

Аннотация

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

Полный текст

Введение

В задачах плазмоники, например при решении нелинейных комплексных дисперсионных уравнений (ДУ), возникает проблема нахождения их комплексных корней. При этом важно построение полной дисперсионной ветви: зависимостей фазовой постоянной (замедления) и постоянной затухания от частоты, например, для поверхностных плазмон-поляритонов (ПП). Мы рассматриваем среды и структуры с диссипацией, поэтому дисперсионные ветви непрерывны и существуют для всех частот, в том числе и в запрещенных для недиссипативных структур зонах. Если в разрешенных зонах при пренебрежимо малой диссипации ПП обычно имеет вид медленной поверхностной волны, а его дисперсия может быть найдена путем решения ДУ в действительной области, то в запрещенной зоне ПП обычно быстрые вытекающие, соответствующие антиповерхностной волне с большим радиационным затуханием. При реальной диссипации ПП классифицируются на втекающие и вытекающие [1‒5], что определяется углом вытекания ϑ=arctank'0x/k'z. Здесь мы считаем, что ПП распространяются вдоль оси z, лежащей на поверхности раздела x = 0 структура – вакуум. Сверху в вакууме k02=k0x2+kz2, k0x=k'0xik''0x, причем мы ищем решение с зависимостью expiωtikzz, kz=k'zik''z. При x<0 может находиться произвольная структура или среда, поддерживающая указанную волну. Волна прямая, если k'zk''z>0 обратная, если k'zk''z<0 [1; 2]. Случай отсутствия диссипации k''z=0 требует для анализа характера волны рассмотрения дисперсии и вычисления групповой скорости. В общем случае определения характера волны следует использовать проинтегрированную по поперечному сечению структуры компоненту вектора Пойнтинга S~z, и тогда критерием служит знак величины S~zk'z [5]. Следует отметить, что для вытекающих волн интеграл S~z может расходиться, и для определения плотности потока мощности его нужно вычислять в определенных границах, а затем переходить к пределу. Для втекания необходимо условие ϑ<0. Определяя Sx, видим, что в этом случае энергия втекает из вакуума в структуру, тогда как при ϑ>0 она вытекает из структуры в вакуум. Границей перехода от медленных к быстрым волнам служит условие «отсечки» n'=k'z/k0=1. Далее будем решать ДУ вида kz=Dk0,kz. Простейшее такое ДУ Ценнека имеет вид kz=k0εω/1+εω и описывает ПП над диэлектрическим полупространством. Оно дает явное решение, поэтому уравнением не является. Для него ПП всегда втекающий, причем для быстрой волны Ценнека при ε'ω>1 втекание может быть весьма слабым, как и для весьма медленного ПП, когда ε'ω1 и диссипация мала: ε''ω/ε'ω1 [4]. В общем случае правая часть ДУ нелинейно зависит от kz даже если возможно описание структуры поверхностным импедансом, поэтому нахождение корней ДУ является актуальной задачей.

Одним из распространенных методов поиска корней общего нелинейного уравнения

fx=0 (1)

является метод простых итерации (МПИ) [6; 7] или последовательных приближений. Согласно ему осуществляются итерации

 xn=φxn1,   n=0,1,2,... (2)

Невязка уравнения (1) совпадает с самой функцией: Δx=fx, а для уравнения (2) невязка есть Δx=xφx. Под невязкой, характеризующей решение, будем понимать квадрат модуля невязки δ=Δx2=Δx2, т. е. квадрат нормы. Функция φx может быть получена из (1) бесконечным числом способов: φx=x+gxfx, где  ‒ gxпроизвольная функция. Обычно для действительного уравнения ее берут непрерывной и знакопостоянной [6]. Если в уравнении (1) функция дифференцируемая p раз, функцию gx также удобно брать таковой. Далее будем предполагать, что все функции дважды, а в некоторых алгоритмах и трижды непрерывно дифференцируемые. МПИ удобен и для комплексного уравнения для определения комплексных корней. Итерационный процесс для (1) и (2) запишем в эквивалентном виде

xn=xn1τnΔn1, (3)

где для уравнения (1) Δn=fxn, а для уравнения (2) Δn=xnφxn. Очевидно, классический МПИ (2) реализуется, если параметр итерации равен единице на каждом шаге: τn=1. Однако выбор параметра итерации на каждом шаге дает дополнительную возможность улучшить сходимость. Это, например, используется в методе минимальных невязок (ММН) [8; 9]: параметр итерации τn выбирается из условия минимума δn на следующем шаге. В случае систем линейных алгебраических уравнений (СЛАУ) или для линейных операторных уравнений A^x=b или x=I^A^x+b в форме (3) xn=xn1τnΔn1, Δn1=A^xn1b имеем [8; 9]:

τn=A^Δn1,Δn1/A^Δn12. (4)

Скалярное произведение понимаем как x,y=x*y. Для сходимости простых итераций (τn1) достаточно, чтобы для одной из норм выполнялось условие I^A^<1. Для сходимости ММН достаточно выполнения условия I^τnA^<1. Считаем матрицу эрмитовой, что всегда можно сделать, умножив СЛАУ на матрицу, эрмитово сопряженную исходной и нормированной: A^=1, разделив СЛАУ на норму. Рассмотрим квадрат спектральной нормы μi. Это наибольшее собственное значение матрицы и по условию μi=1. Тогда итерации выполняются с параметром (4), если 0<τn<1. При выходе из этой области параметр следует корректировать. При τn>1 следует брать τn=1δ, 0<δ<1. При очень малом параметре алгоритм почти зацикливается, поэтому следует брать τn=δ, не допуская слишком малых значений. Алгоритм метода скорейшего спуска (МСС) аналогичен, но в нем минимизируется энергетическая норма невязки, а параметр итерации приобретает известный вид τn=n12/A^Δn1,Δn1. ММН и МСС могут не сойтись. Даже в случае одного линейного уравнения итерационный процесс может разойтись, а может сходиться достаточно медленно. Для одного линейного уравнения ММН и МСС дают решение за один шаг в силу скалярной природы задачи: τ1=A1. Это также следует из условия Δ12=0. В многомерном случае Δn20, и из наложения условия Δn2=0 не следует, вообще говоря, сходящийся алгоритм. Получающееся квадратное уравнение для τn может иметь комплексные значения для действительной СЛАУ. Использование матричного параметра итерации τ^n сопряжено с трудностями его определения: такой оптимальный параметр должен быть взят в виде τ^n=A^1, т. е. приводит к той же самой задаче, для решения которой он и выбирается. Матричный параметр целесообразен в специальных случаях, например, когда матрица A^ имеет большие диагональные элементы. Однако для поиска корня нелинейного уравнения условие Δn2=0, весьма продуктивно.

  1. Методы коррекции параметра итерации

При решении нелинейного скалярного уравнения соотношение Δn2=0 может быть наложено, поскольку принципиально может быть достигнуто. Цель работы – модификация МПИ с целью получения сходящихся алгоритмов для решения нелинейных уравнений. В работе [9] рассмотрен ММН в случае расположения начальной точки x близко к корню x¯. Он основан на разложениях Тейлора типа fxf'x¯xx¯+f''x¯xx¯2/2, в частном случае – на линеаризации, т. е. на удержании в разложениях такого типа линейного члена. Рассматривая уравнение (2) в форме (3), имеем δn=xnxn1+τnΔn12. Параметр τn можно искать из условия, что δn имеет минимум, т. е. из условия δn/τn*=0 [8]. Поскольку Δn1*0 (иначе итерации уже бы сошлись на шаге n1), это эквивалентно условию Δn=0. Из него имеем τn=xn1xn/Δn1. Поскольку Δn=0, это означает, что xn=x¯=φx¯. С другой стороны, φxn1=φx¯+xn1x¯x¯+φ'x¯xn1x¯, поэтому имеем соотношение

Δn=x¯xn1+τnxn1x¯φ'x¯xn1x¯,

из которого следует τn=1φ'x¯1. Итерационный процесс приобретает вид

xn=xn111φ'xn1++1φ'xn1φxn1.

Поскольку корень  неизвестен, процесс реализуем заменой xnxn1 в форме

xn=φ~xn1=1τnxn1+τnφxn1, (5)

полагая τn=1φ'xn11 Итерации (5) имеют форму (2) с переопределением функции φ на φ~x. Рассмотрим, к чему он приводит. Известно, что процесс (2) сходится, если φ'xq<1 в области сходимости (принцип сжимающих отображений). Однако если это условие не выполнено, процесс (2) может разойтись. Это требует поиска и подбора вида функций φx. Удобно строить такие процессы, которые сходятся из любого начального приближения. Конечно, теоретическая сходимость не означает возможности вычисления. Скажем, вполне можно построить процесс, сходящийся за несколько итераций при движении из начальной точки вблизи корня, и сходящийся за 106 или за существенно большее число итераций при движении из удаленной точки. Такие точки нецелесообразно рассматривать и всегда желательно определять области, из которых скорость сходимости высокая. Если рассмотреть форму (5), то для выполнения условия сжимающих отображений достаточно выполнить 1τn1φ'xn1<1/α, где параметр α>1. Выбирая этот параметр на каждом шаге, можно находить параметр итерации и строить сходящийся итерационный процесс. Удобнее, однако, сразу взять α, т. е. положить τn=1φ'xn11. В этом случае производная правой части (5) в точке xn1 точно равна нулю, а процесс сходится, даже если φ'xn1>1. Такие образом, на каждом шаге процесс (5) подбирает оптимальную функцию φ~x в правой части так, что φ'~xn1=0.

Рассмотрим ряд частных случаев. Если φ'xn1=0, то τn=1, в окрестности этой точки реализуется хорошо сходящийся процесс (2). Если φ'xn11 (точка вблизи полюса производной), то τn0, поэтому xnxn1. Это медленно сходящийся процесс, и для его убыстрения целесообразно сделать одну-две простые итерации (2), а затем перейти к процессу (5) в новой точке, т. е. отойти от полюса производной. Если φ'xn1=xn1, то xn1xn1=xn1xn1+φxn1. Вблизи корня функция отличается от аргумента на малую величину: φxn1=xn1+δ. Тогда xn=xn1+δ/1xn1. Если xn1 существенно меньше единицы или величина xn1 отрицательная, то сходимость хорошая. Если φ'xn11 возможна разболтка. Если φ'xn1<1 то для ее преодоления надо сделать 2‒3 итерации (2), а если φ'xn1>1 то 2‒3 итерации (5) без пересчета параметра итерации, а затем продолжить процесс (5). Если 1φ'xn1δ, где величина δ порядка точности представления чисел в ЭВМ (погрешность округления), это означает, что точка вблизи корня, в котором φ'x¯=1 При достижении такой точки целесообразно фиксировать τn и выполнять метод простой итерации с τn1 Вообще, если φ'xn1=1 это означает, что xn1=φxn1, т. е. xn1 ‒ это корень. Поэтому рассматриваем случай φ'xn11 Это может быть, только если φ'x¯=1 и процесс подошел к корню. Фактически продолжение итераций (5) означает использование медленно сходящегося процесса (2). В этом случае следует сделать несколько итераций (5) без пересчета параметров, что обычно достаточно для получения результата. При приближенном вычислении производной из-за разболтки сходимость может не достигаться. В этом случае удобно вообще зафиксировать параметр итерации, пока производная слабо отличается от единицы. Если φ'x¯=1 для достижения быстрой сходимости функцию целесообразно заменить на новую φ~x=x+gxxφx. На эту функцию наложим условие φ'~xn1=0 что делает процесс быстро сходящимся. Функцию gx можно взять с точностью до любого множителя, поэтому gxn1 задаем произвольно. Полагаем φ'~x=1+g'xxφx+gx1φ'x=axxn1.

Справа вообще можно написать полином, но это усложняет задачу. Имеем 1+g'xn1 xn1φxn1=0, φ~x=axxn12/2+φ~xn1.

Параметр a возьмем из условия, что сходимость в точке начального приближения хорошая, например, берем φ'~x0=1/2. Таким образом, необходимо интегрировать уравнение

g'x=axxn11gx1φ'xxφx

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

gxn1=φ~xn1xn1xn1φxn1.

Интегрирование желательно выполнить аналитически, например методом рядов. Если это возможно, к новой функции применимы оба процесса (2) и (5). Начальное условие для gx произвольное и может использоваться для упрощения решения, как и параметр a.

Поскольку при итерации (5) xn не попадает точно в корень x¯, итерации продолжаются, однако погрешности xnxn1 не возрастают, а обычно уменьшаются. Численное вычисление производных может приводить к разболтке, т. е. к немонотонному изменению xnxn1. Вблизи корня процесс имеет второй порядок и сходится квадратично, т. е. подобно методу Ньютона [6]. Однако метод Ньютона сходится не всегда. Преимущества процесса (5) в том, что по построению φ'~xn1=0, n=0,1,2,...,  поэтому он обладает бесконечной областью гарантированной сходимости. Однако при его реализации и численном вычислении производных возможно зацикливание. Его определяем условиями, когда невязка xnxn1 мала или изменяется немонотонно, а невязка уравнения Δn еще не мала для выхода из процесса. Тогда при численном вычислении производных для избегания зацикливания нужно после итерации (5) провести 2‒3 простые итерации (2), если для них φ'<1, или по формуле (5) без пересчета параметра итерации, если φ'>1. С использованием полученного параметра τn и функции fx=xφx метод Ньютона в форме xn=xn1xn1φxn1/1φ'xn1 приобретает вид (5), т. е. сходится, если φ''xn1<1φ'xn12/Δn1. Если вторая производная в корне ограничена, то сходимость имеет место всегда при попадании в малую окрестность корня. Быстрая сходимость в окрестности корня имеет место, если φ''x¯=0.

 

Таблица. Сходимость (ход итераций) процессов для вычисления x=2

Table. Сходимость (ход итераций) процессов для вычисления x=2

n

x=x/2+1/x

x=4x+2/x1 и алгоритм (5)

Алгоритм

(10)

x=x/2+1/x

x=4x+2/x1

Алгоритм

(10)

0

10,0

10,0

10,0

0,1

0,1

0,1

1

5,100000

0,392157

6,733333

10,050000

0,199005

6,733333

2

2,746078

0,728311

4,587899

5,124503

0,390282

4,587899

3

1,737195

1,151281

3,203909

2,757392

0,725323

3,203909

4

1,444238

1,3841813

2,344018

1,741358

1,148529

2,344018

5

1,414526

1,413902

1,847091

1,444943

1,384137

1,847091

6

1,414214

1,414214

1,592322

1,414540

1,413887

1,592322

 

Рассмотрим в качестве примера классическую задачу извлечения квадратного корня x=a. Как известно, алгоритм x=a/x всегда расходится из любой точки, а алгоритм x=x+a/x/2 сходится всегда из любого начального приближения [6]. Причина в том, что в первом случае φ'x=a/x2, φ'x1 при малых x2, и в точке корня φ'x¯=1, а во втором случае φ'x=1a/x2/2, и в точке корня φ'x¯=0. При этом всегда имеет место попадание в малую окрестность корня. Если в методе Ньютона взять fx=x2a, получим сходящийся алгоритм x=x+a/x/2 с τn=1/2. Если в методе Ньютона взять fx=x2a, получим расходящийся алгоритм x=a/x. Вообще для функции fx=x2aα сходимость метода Ньютона будет, если 1α1<1. Следуя вышеизложенному, формируем процесс x=φ~x=1τx+aτ/x, в котором τ=1+a/x21. Как видно, φ~x=4x+a/x1. Для этой функции φ'~a=0. Желательно использовать алгоритм (5), для которого φ'~xn=0, n = 1, 2, … Он эквивалентен алгоритму (2) с приведенной функцией φ~x. Однако в общем случае, когда производная неизвестна и вычисляется численно, получить функцию φ~x нельзя. В случае извлечения корня такая формула приобретает вид x=2x+a/x/3. Для нее, как и для процессов x=x/2+1/x, x=4/x+2/x и (5), в таблице приведены сходимости итераций для двух начальных приближений x0=10 и x0=0,1 при a=2. Еще один процесс с удержанием производной второго порядка следует из равенства xn=φxn1τnΔn1 при условии, что величина ξn=τnΔn1 мала по модулю. Получаем

φxn1ξnφxn1φ'xn1ξn+φ''xn1ξn2/2,

что дает

ξn=ξxn11φ'xn1/φ''xn1++1φ'xn1/φ''xn12+2Δn1/φ''xn1.

  1. Решение дисперсионных уравнений для поверхностных плазмонов

Поверхностные ПП в настоящее время широко исследуются в нанооптике, фотонике и наноплазмонике [1‒5; 10‒19]. Особенность их моделирования в том, что необходим поиск комплексных корней нелинейных трансцендентных дисперсионных уравнений. Рассмотрим решение комплексного уравнения

n=1z2n. (6)

Это уравнение описывает комплексный коэффициент замедления ПП вдоль поверхности образца с комплексным импедансом zn. Поверхностный импеданс характеризует входное сопротивление образца. Металл описывается комплексной зависящей от круговой частоты  диэлектрической проницаемостью (ДП) εω=εLωp2/ω2iωωc, где ωp, ωc ‒ константы, а εL на низких частотах также является константой. На высоких частотах мы опишем ее двумя термами дисперсии Лоренца, подбирая параметры так, чтобы удовлетворить экспериментальным данным для металлов в оптическом диапазоне. Далее обозначаем ε=ε'iε''. Уравнение (6) имеет точное решение n=ε/ε+1 для любого однородного образца с ДП e, занимающего полупространство (J. Zenneck, 1907). Это точное решение получается сшиванием электромагнитных полей или приравниванием входных импедансов, например металла zn=εn2/ε и характеристического импеданса вакуума z0n=1n2  для E-волны. Однако если образец неоднородный (слоистый или ДП ε зависит от нормальной координаты), входной импеданс zn является сложной функцией, но уравнение (6) все равно справедливо. Решения таких уравнений рассмотрены в ряде работ (см. [1‒4; 14]). Кроме того, в ряде работ показано совместное итерационное решение интегральных уравнений и соответствующих им функционалов, получающиеся из линейной теории рассеяния или линейной электродинамики. Искомый параметр в такие функционалы входит нелинейно, в том числе и через аргументы трансцендентных функций. В нелинейной теории рассеяния итерационные методы, как правило, являются наиболее эффективными для получения решений. Часто в теории рассеяния возникают интегральные уравнения типа Липпмана – Швингера в форме объемных (трехмерных) уравнений Фредгольма второго рода, для которых удобно организовать итерационный процесс типа (2), в котором φ^x ‒ интегральный или интегродифференциальный оператор, а x ‒ функция трехмерных координат и времени. Для таких задач также применимы изложенные итерационные алгоритмы.

Для задач о плазмонах есть области, где модуль производной правой части уравнений типа (6) больше единицы и метод прямой итерации расходится. Решение n=ε/ε+1 для металла (серебро) приведено на рис. 1, кривая 1. Однако если решать уравнение (6) МПИ, то в области ε'<1 он сходится, в области 1<ε'<1 расходится, а в области 1<ε' снова имеет место сходимость. Действительно, φn=n2+ε2ε/ε, φ'n=n/εn2+ε2ε. Будем считать потери очень малыми, т. е. ε''ε'. В области малых частот n1, ε21, и φ'n0. В области плазмонного резонанса ε'1 имеем n1i/2ε'', n1 и φ'n1. Если ε0, то n0, а φ'n1. Если ε=1, то n=1 и φ'1=1. Нетрудно видеть, что φ'n=1/ε2=1n22/n4. Поэтому сходимость возникает при ε>1, когда 1/2<n<1 и φ'n=1/ε2<1. Построим процесс, который сходится во всей области комплексных значений ДП. Если взять φ'n=1/ε2, процесс типа (2) определен функцией

φ~n=ε2n2/ε2+11/εn/ε21. (7)

Однако такой процесс хорошо сходится только в окрестности корня: именно там φ'~n=0. Процесс с индексом k, сходящийся из любого начального приближения, имеет вид

nk=nk11τk+τknk12/ε2+11/ε. (8)

Для него nk/nk1=0, τk/nk1=0, где

τk=ε2nk2/ε2+11/εε2nk2/ε2+11/ε1.

Реально необходимо строить частотную дисперсию ПП. Полагая для низких частот n0=1, получаем быстро сходящийся к решению n1 процесс. Это решение на низких частотах определяет быстрый (движущийся почти со скоростью света) ПП. Плавно повышая частоту и используя для начального значения полученное решение на предыдущей частоте, получаем быстро сходящиеся процессы для всех частот. Преимущество таких вычислений в том, что всегда начальное приближение близко к корню, поэтому скорость сходимости высокая, и можно также использовать процесс (7). Результаты расчета приводят к точкам 2 на рис. 1, точно ложащимся на дисперсионную ветвь (кривая 1, рис. 1).

 

Рис. 1. Дисперсия замедления ПП на поверхности серебра (ωp=1,351016, ωc=2,71013, εL=9,0) (кривая 1), на поверхности серебра, покрытой диэлектрической пленкой толщины 250 нм с ДП 4i0,004 (3), на поверхности серебра, покрытой диэлектрической пленкой и пленкой серебра толщиной 10 нм (4). Точки 2 показывают сходимость итераций (5) для ПП на серебре

Fig. 1. Dispersion of PP deceleration on the silver surface (ωp=1,351016, ωc=2,71013, εL=9,0) (curve 1), on the silver surface covered with a dielectric film 250 nm thick with DP 4i0,004 (3), on the surface of silver coated with a dielectric film and a silver film 10 nm thick (4). Points 2 show the convergence of iterations (5) for PP on silver

 

Рассмотрим более сложную подобную задачу. Пусть на металле лежит слой диэлектрика с ДП εd толщины d. Тогда входной импеданс приобретает вид

zn=zdz0n+izdtank0dεdn2zd+iz0ntank0dεdn2. (9)

Опять решаем уравнение (6). Результат решения приведен на рис. 1, кривая 3. Еще один пример на рис. 1 (кривая 4) соответствует случаю, когда еще и на слое диэлектрика лежит слой металла толщины dm. Импеданс строится рекуррентно по формуле (9) подстановкой в нее z0nzn и zdzm, εdε, ddm. Можно получить метод более высокого порядка, оставляя в разложении правой части в точке корня три члена φx=φx¯+φ'x¯xx¯+φ''x¯xx¯2/2 и полагая δn/τn*=0, что также ведет к условию Δn=0. Окончательно имеем

x=x2φ'xφ'''xφ''2xφxφ''2x2φ'xφ'''x. (10)

 

Рис. 2. Дисперсия потерь ПП на поверхности серебра (1), на серебре, покрытом диэлектрической пленкой (3), на серебре, покрытом диэлектрической пленкой и пленкой серебра толщиной 10 нм (4). Точки 2 показывают сходимость итераций (5) для ПП на серебре. Структура соответствует рис. 1

Fig. 2. Dispersion of PP losses on the surface of silver (1), on silver coated with a dielectric film (3), on silver coated with a dielectric film and a silver film with a thickness of 10 nm (4). Points 2 show the convergence of iterations (5) for PP on silver. The structure corresponds to Fig. 1

 

Рис. 3. Дисперсия в слое асимметричного ГММ (вкладка) толщины 50 нм при разных углах a: p/3, p/4 и p/6. Знаки «+» и «‒» в y соответствуют графикам (а) и (б)

Fig. 3. Dispersion in an asymmetric HMM layer (inset) 50 nm thick at different angles a: p/3, p/4, and p/6. The signs «+» and «‒» in y correspond to graphs (a) and (b)

 

Рассмотрим пример более сложной структуры в виде слоя гиперболического метаматериала ГММ второго типа [5] (см. вставка на рис. 3, а) толщины d, образованного тонкими серебряными пленками, разделенными диэлектрическими слоями SiO2. Для упрощения гомогенизации мы считаем слои металла и диэлектрика одинаковой наноразмерной толщины. Тогда гомогенизация справедлива вплоть до оптических частот, и слой, имеющий ось анизотропии вдоль оси z, характеризуется диагональным тензором ДП εxx=εyy=εm+εd/2, εzz=2εmεd/εm+εd.  Здесь мы обозначили ДП металла εm и диэлектрика εd. Для ДП металла (серебра) мы используем модель Друде – Лоренца с введением двух Лоренцевых термов для описания сложного поведения ε'm вблизи перехода через нуль [20]. Если повернуть ось анизотропии на угол α относительно оси y, получим тензор ДП вида

ε^=cos2αεxx+sin2αεzz00εxxsinαcosαεzzεxx0       sinαcosαεzzεxx0cos2αεzz+sin2αεxx. (11)

Он соответствует случаю асимметричного ГММ, когда ось анизотропии направлена под углом к границе раздела. Еще более сложным является случай, когда волна распространяется вдоль поверхности в произвольном направлении, т. е. образует некий угол q с осью z. Поворачивая систему координат вдоль оси x так, что новая ось z совпала с направлением распространения, получим в новой системе координат тензор ДП вида

ε~=ε~xxε~xyε~xzε~xyε~yyε~yzε~xzε~yzε~zz, (12)

с компонентами:

ε~xx=cos2αεxx+sin2αεzz,

ε~yy=cos2θεxx+sin2θ××cos2αεzz+sin2αεxx,

ε~zz=sin2αεxx+cos2θ××cos2αεzz+sin2αεxx,

ε~xy=sinθsin2αεzzεxx/2,

ε~xz=cosθsinαcosαεzzεxx,

ε~yz=sin2θ××εxxcos2αεzz+sin2αεxx/2.

Мы не приводим довольно сложно получающиеся ДУ. Наиболее простой случай здесь соответствует полупространству с тензором ДП (11), когда волны распадаются на E и H ПП. Также имеется два типа волн для случая конечного слоя (рис. 3). В этом случае ДУ имеет вид

kz±=±k0ε~xxfk0,kz/Δ1fk0,kz/Δ11/2, (13)

где Δ=ε^xxε^zzε^xz2, ψ=±dk02ε^xxkz2Δ/ε^xx, а входящая в (13) функция представляется в виде

fk0,kz==ρρ0expiψρ+ρ0expiψρρ0expiψ+ρ+ρ0expiψ2.

В случае тензора ДП (12) волны вдоль поверхности слоя гибридные. ДУ может быть получено методом сшивания и на основе матриц передачи (переноса)  – так называемых матриц Берремана. Для рис. 3 также существуют области прямых и обратных волн на дисперсионных кривых. В силу корневых многозначностей в общем случае возможны четыре дисперсионные ветви. Возможно рассмотрение ПП вдоль поверхностей ГММ, построенных на основе графен-диэлектрических структур [5]. ГММ с полупроводниковыми и графеновыми слоями хорошо оптически управляются. Также возможно создание гиротропных структур, управляемых магнитным полем, но получение ДУ в этом случае является весьма сложной задачей.

Заключение

Показано, что алгоритм МПИ с коррекцией сходимости путем подбора правой части, имеющей нулевую производную на каждом шаге, совпадает с ММН первого приближения (при линеаризации) и имеет то преимущество, что всегда сходится (по крайней мере, невязка не возрастает). Такой выбор τn приводит к равенству нулю линеаризованной невязки. Вблизи корня сходимость квадратичная. Однако в ряде случаев возможно сильное замедление сходимости. В этих случаях можно предусмотреть специальные методы коррекции шага. Если корней несколько, алгоритм их поиска может состоять из движений из сильно различных начальных приближений.

Рассмотренные адаптивные методы удобны для задач плазмоники и позволяют определять целиком всю комплексную ветвь и даже несколько ветвей (если они есть) дисперсионного уравнения при плавном изменении частоты (или иного параметра), т. е. решать проблему поиска и классификации комплексных корней. На рис. 2 приведены результаты для потерь n''ω, соответствующие рис. 1. Приведенные численные результаты, рис. 1–3 показывают возможность существования обратных ПП, для которых n'ωn''ω<0. Если дисперсионная кривая 4, рис. 1 строится из условия Renω>0, то обратным ПП соответствуют «отрицательные» потери (рис. 2, кривая 4).

×

Об авторах

Михаил Владимирович Давидович

Саратовский национальный исследовательский государственный университет имени Н.Г. Чернышевского

Email: davidovichmv@info.sgu.ru

доктор физико-математических наук, профессор кафедры радиотехники и электродинамики 

Россия, Саратов

Александр Константинович Кобец

ООО «НПФ “Этна плюс”»

Email: kobetzak@info.sgu.ru

директор по качеству 

Россия, Саратов

Кирилл Александрович Саяпин

Саратовский национальный исследовательский государственный университет имени Н.Г. Чернышевского

Автор, ответственный за переписку.
Email: sayapin_kirill@mail.ru

аспирант кафедры радиотехники и электродинамики Института физики

Россия, Саратов

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

  1. Давидович М.В. Плазмоны в многослойных плоскослоистых структурах // Квантовая электроника. 2017. Т. 47, № 6. С. 567‒579. URL: http://mi.mathnet.ru/qe16620
  2. Давидович М.В. Втекающие и вытекающие несобственные моды – анализ диссипативных дисперсионных уравнений и волна Ценнека. Саратов: Сарат. ун-т, 2014. 104 с.
  3. Давидович М.В. Максимальное замедление и отрицательная дисперсия плазмонов вдоль металлического слоя // ПЖТФ. 2017. Т. 43, вып. 22. С. 55‒62. DOI: https://doi.org/10.21883/PJTF.2017.22.45261.16629
  4. Давидович М.В. Об условии перехода быстрой поверхностной волны в медленную // Радиотехника и электроника. 2018. Т. 63, № 6. С. 499–506. URL: https://www.elibrary.ru/item.asp?id=35104331
  5. Давидович М.В. Гиперболические метаматериалы: получение, свойства, применения, перспективы // УФН. 2019. Т. 189, № 12. С. 1250‒1284. DOI: https://doi.org/10.3367/UFNr.2019.08.038643
  6. Калиткин Н.Н. Численные методы. М.: Наука, 1978. 512 с.
  7. Chapra S.C., Canale R.P. Numerical Methods for Engineers. New York: McGraw-Hill Education, 2014. 970 p.
  8. Самохин А.Б. Интегральные уравнения и итерационные методы в электромагнитном рассеянии. М.: Радио и связь, 1998. 160 с.
  9. Давидович М.В. Итерационные методы решения задач электродинамики. Саратов: Сарат. ун-т, 2014. 240 с.
  10. Long-range plasmonic waveguides with hyperbolic cladding / V.E. Babicheva [et al.] // Opt. Express. 2015. Vol. 23, no. 24. P. 31109–31119. DOI: https://doi.org/10.1364/OE.23.031109
  11. Lyashko E.I., Maimistov A.I. Guided waves in asymmetric hyperbolic slab waveguides: the TM mode case // J. Opt. Soc. Am. B. 2016. Vol. 33, no. 11. P. 2320–2330. DOI: https://doi.org/10.1364/JOSAB.33.002320
  12. Chern R.-L., Yu Y.-Z. Chiral surface waves on hyperbolic-gyromagnetic metamaterials // Opt. Express. 2017. Vol. 25, no. 10. P. 11801–11812. DOI: https://doi.org/10.1364/OE.25.011801
  13. Chiral surface waves supported by biaxial hyperbolic metamaterials / W.-L. Gao [et al.] // Light Sci. Appl. 2015. Vol. 4, no. 9. P. e328-1–5. DOI: https://doi.org/10.1038/lsa.2015.101
  14. Popov V., Lavrinenko A.V., Novitsky A. Surface waves on multilayer hyperbolic metamaterials: Operator approach to effective medium approximation // Phys. Rev. B. 2018. Vol. 97, no. 12. P. 125428-1–10. DOI: https://doi.org/10.1103/PhysRevB.97.125428
  15. Engineered surface Bloch waves in graphene-based hyperbolic metamaterials / Y. Xiang [et al.] // Opt. Express. 2014. Vol. 22, no. 3. P. 3054–3062. DOI: https://doi.org/10.1364/OE.22.003054
  16. Tunable terahertz amplification based on photoexcited active graphene hyperbolic metamaterials [Invited] / T. Guo [et al.] // Opt. Materials Express. 2014. Vol. 8, no. 12. P. 3941–3952. DOI: https://doi.org/10.1364/OME.8.003941
  17. Grig T., Hess O. Tunable surface waves at the interface separating different graphene-dielectric composite hyperbolic metamaterials // Opt. Express. 2017. Vol. 25, no. 10. P. 11466–11476. DOI: https://doi.org/10.1364/OE.25.011466
  18. Tunable bulk polaritons of graphene-based hyperbolic metamaterials / L. Zhang [et al.] // Opt. Express. 2014. Vol. 22, no. 11. P. 14022–14030. DOI: https://doi.org/10.1364/OE.22.014022
  19. Давидович М.В. Плазмон-поляритоны Дьяконова вдоль гиперболического метаматериала // Компьютерная оптика. 2021. Т. 45, № 1. С. 48–57. DOI: https://doi.org/10.18287/2412-6179-CO-673
  20. Johnson P.B., Christy R.W. Optical constants of the noble metals // Phys. Rev. B. 1972. Vol. 6, no. 12. P. 4370–4379. DOI: https://doi.org/10.1103/PhysRevB.6.4370

Дополнительные файлы

Доп. файлы
Действие
1. JATS XML
2. Рис. 1. Дисперсия замедления ПП на поверхности серебра (кривая 1), на поверхности серебра, покрытой диэлектрической пленкой толщины 250 нм с ДП (3), на поверхности серебра, покрытой диэлектрической пленкой и пленкой серебра толщиной 10 нм (4). Точки 2 показывают сходимость итераций (5) для ПП на серебре

Скачать (131KB)
3. Рис. 2. Дисперсия потерь ПП на поверхности серебра (1), на серебре, покрытом диэлектрической пленкой (3), на серебре, покрытом диэлектрической пленкой и пленкой серебра толщиной 10 нм (4). Точки 2 показывают сходимость итераций (5) для ПП на серебре. Структура соответствует рис. 1

Скачать (119KB)
4. Рис. 3. Дисперсия в слое асимметричного ГММ (вкладка) толщины 50 нм при разных углах a: p/3, p/4 и p/6. Знаки «+» и «‒» в y соответствуют графикам (а) и (б)

Скачать (389KB)

© Давидович М., Кобец А., Саяпин К., 2021

Creative Commons License
Эта статья доступна по лицензии Creative Commons Attribution 4.0 International License.

СМИ зарегистрировано Федеральной службой по надзору в сфере связи, информационных технологий и массовых коммуникаций (Роскомнадзор).
Регистрационный номер и дата принятия решения о регистрации СМИ: серия ФС 77 - 68199 от 27.12.2016.

Данный сайт использует cookie-файлы

Продолжая использовать наш сайт, вы даете согласие на обработку файлов cookie, которые обеспечивают правильную работу сайта.

О куки-файлах