Объединенный метод конечных элементов и граничных элементов для анализа дифракции света на дифракционных решетках

ОБЪЕДИНЕННЫИ МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ И ГРАНИЧНЫХ ЭЛЕМЕНТОВ ДЛЯ АНАЛИЗА ДИФРАКЦИИ СВЕТА НА ДИФРАКЦИОННЫХ РЕШЕТКАХ

Д.В.Нестеренко1, В.В. Котляр1,2 1 Институт систем обработки изображений РАН,

2 Самарский государственный аэрокосмический университет имени академика С.П. Королева

Аннотация

Рассматривается применение объединенного метода на основе конечных элементов и граничных элементов для анализа задачи дифракции электромагнитной волны на цилиндрическом периодическом объекте, период которого сравним с длиной волны. Результаты моделирования диэлектрических дифракционных решеток с различным периодом были сравнены с аналогичными результатами, полученными известным методом связанных волн.

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

Введение

Теория рассеяния на периодических структурах, обычно называемых дифракционными решетками, имеет много применений в микрооптике, например, электромагнитные и оптические средства связи, средства визуализации, определение свойств объектов и поверхностей, электронные и оптические компоненты, фотонные кристаллы, дифракционные решетки

[1]. Для моделирования дифракции света на дифракционных решетках были разработаны численные методы, среди которых известны дифференциальные и интегральные методы, методы, основанные на распространении Рэлеевских и собственных мод, вариационные и конечно-разностные методы: метод связанных волн

[2], С-метод [3], методы конечных элементов [4, 5, 6], интегральные методы [7], разностные (finite difference-time domain, FDTD) методы [8, 9].

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

В качестве частного случая вариационных методов может быть рассмотрен метод конечных элементов (МКЭ), примененный к эллиптическому уравнению Гельмгольца в области расчета. Он включает выбор схемы дискретизации, построение и минимизацию функционального соотношения. Полученное соотношение преобразуется в систему линейных уравнений, являющуюся неполной без применения граничных условий.

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

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

В настоящей работе описывается формулировка объединенного метода для задач рассеяния света на периодических объектах на основе метода конечных элементов и метода граничных элементов (ПМКЭ-ГЭ). С помощью разработанного метода ПМКЭ-ГЭ и RCWA метода [10] было проведено сравнительное моделирование дифракции света на одномерной диэлектрической дифракционной решетке. Сравнение результатов моделирования приведено для ТЕ- и ТМ-поляризованных волн.

Дифракция на периодических структурах

Рассмотрим дифракцию плоской волны с волновым вектором k = k (sin (в), - cos (в) ,0), k = k04s

на диэлектрической периодической структуре с периодом d, k0 - волновое число волны в свободном пространстве

k 2Л

0 = V

где ^ - длина волны в свободном пространстве, s - диэлектрическая проницаемость среды.

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

как сумма рассеянного и падающего полей. Цель задачи - определить амплитуду и фазу отраженных, прошедших и затухающих порядков.

Геометрия задачи приведена на рис. 1, где Rn и Тп - коэффициенты отражения и пропускания дифракционных порядков. Для приведенной геометрии задачи можно выделить три зоны с различной диэлектрической проницаемостью: область над структурой при у > а (^), где а - максимальная высота структуры с постоянной диэлектрической проницаемостью е= еь область структуры 0 < у < а с диэлектрической проницаемостью е= е(х, у) и область У < 0 (Т3) с постоянной диэлектрической проницаемостью е= е3.

Задача дифракции плоской волны на одномерной периодической структуре сводится к рассмотрению двух независимых задач: задачи дифракции плоской волны с ТЕ-поляризацией (Е2 Ф 0, Нг = 0) и задачи дифракции плоской волны с ТМ-поляризацией (Н Ф 0, Е = 0) [11].

Полное поле иО(х, у) в области О (0 < х < d, 0 < у < а) должно удовлетворять уравнению [12]:

V'

1

-^п

+^q

_ р

где /а = jk0Z0 Jz, р(х, у)=Мг, q(x, у)=—г для ТЕ - поляГ (л °У1

• г, р(х, у)=єг, q(х, у)=/иг

и /а =- (л

Vx — 1

V—

для ТМ - поляризации. Константы /иг и ег представляют собой отношение магнитной и диэлектрической проницаемостей среды к аналогичным показателям свободного пространства, т.е. и^иио и ег=е/е0, =^и0 / е0 - импеданс свободного про-

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

Для решения уравнения (1) область расчета О должна быть ограничена введением искусственной границы Г = Г + Г2 + Г3 + Г4 (см. рис. 1). Г и Г3 -фиктивные границы, бесконечно простирающиеся параллельно оси х по координатам у = а и у = 0. Соответственно, для единственного решения задачи на данной искусственной границе должны быть введены граничные условия.

Так как пространство в зонах и Т3 однородное, то поле в них может быть определено в терминах граничных интегралов с соответствующей функцией Грина. Полное поле и^(х, у) в этих зонах должно удовлетворять уравнению:

V

- Vuр (#) Р

+ ^ит (#) = /р ,

где /р = jk0Z0л ,Р = Мт, q = —т

для ТЕ-поляризации и /р = -

Vy.

• г, р = єт,

q = м т для ТМ - поляризации. Л - вектор плотности электрического тока источника в области Р, область Р = ¥ь Р3.

Рис. 1. Геометрия задачи дифракции на периодической диэлектрической структуре

Описание метода расчета Метод решения уравнения Галеркина (1) основан на решении соотношений вида:

Я

где у - произвольная функция из области определения уравнения (1).

Используя первую формулу Грина:

Л Р^а О = | Р ^^1 - Л VPVQd О,

О г ап О

для функций Р и Q, где О - область плоскости х, у;

Г - ее граница, обходимая против часовой стрелки;

ап

к кривой Г, получим:

производная в направлении внешней нормали

я[ р уи^ (X. у№у- qk 2иа (х, у)г|/ -

-І

г .

^ ёиа (X, у) ёа

ёГ = ([ /аудП.

Систему базисных функций для  обозначим

N ,Ыу 0

{а>Пк,і(х, у)} =0у и систему базисных функций для Г

обозначим {а>гт(х,у)}где Ш, ^ - число узлов сеточного покрытия прямоугольной области О по

оси х и у соответственно, М - число узлов сеточного покрытия границы Г.

Подставляя в соотношение (4) вместо произвольной функции у систему базисных функций для метода Галеркина, можно записать систему линейных уравнений:

Аи + Bv = С/

где и = (щ, ..., и ЬГхЫу )т - вектор, составленный из коэффициентов (к)+1 ; разложения:

\Xу) = 2 uk,іфи(х,у)-

k ,1=0

Вектор f = (/1, ..., / МхКу) - вектор, составленный из коэффициентов разложения:

х ' у

/(х, у) = 2 /кіію1і(х, у).

Хотя равенства (6) и (7) действительны для всех точек (х, у) в области , необходима отдельная обработка величин поля и его частных производных на границе Г от значений во внутренней области. Разложение, аналогичное (6) и (7), для поля и его частных производных на границе имеет вид:

и Г( х у) = 2 иХ(х y),

т=1

М

VГ(X, у) = 2 VnPYm ( X, у) =

т=1

М

/Г(X, у) = 2 /т<Г (X, У),

т=1

где (х, у) е Г, V = (VI, ..., V") - вектор, составленный из коэффициентов разложения vk.

Элементы матрицы А вычисляются по формулам:

aNy (К) +1, Ny (/) + j =

=я

1

Р

+

- +

дх дх

да%,

ду ду

-ko2 q(x, у)ю° (х, у)5 (х, у)) а О,

К, I = [1, N1; I,] = [1, Ny], где ^ - область разбиения области ^, включающая узлы сети К и].

Элементы матрицы В вычисляются по формулам:

Ьт,, =- <? ^5 ё1-.

Гт 5

т, 5 = [1, М],

где Гт, 5 - линейная область границы Г, включающая узлы границы т и 5.

Элементы матрицы С вычисляются по формулам:

С

= Ца?,1

К, I = [1, N1; I,} = [1, Ny], где j - область разбиения области ^, включающая узлы сети К иj.

В качестве кусочно-линейного базиса была определена функция вида:

хк - х у - у

(х> у) =

1^^^— h

1-х-х h ’

у - у

h

если X, у єС^ если X, у Є hi2

1+

h ’

если х

13

,

х—х + у—у, если х, у єО?,4

hh %14

1+, h

уі - у

1-

h

если х, у є0^5

если х, у Є^

где h - длина сегмента сети покрытия.

Элементы (а1к\) матрицы А, элементы (Ьт,5) матрицы В и элементы (с К";) матрицы С вычисляются

из уравнений (11), (12) и (13) соответственно. Тогда систему уравнений (5) можно представить в виде:

IА0,0] [Аг,а ] 0 ] иа

_[А0,г] [Аг,г ] [в]. иг _ V

[С^ГСгЛ С,г ІСг,г]

г

Система уравнений (15) не имеет единственного решения, так как она состоит из N равенств с Ы+М неизвестными: N = ЫхЫу - общее число узловых величин поля ик1(х, у) в области  и М производных по нормали на граничных узлах ук>1(х, у).

Определим граничные условия для поля и его производных на границах Гі и Г3, удовлетворяющие уравнению (2).

Для построения граничного интегрального уравнения на основе уравнения (2) для поля и его нормальной производной введем функцию Грина и , удовлетворяющую условиям излучения Зоммер-фельда и являющуюся фундаментальным решением уравнения Гельмгольца:

V2и * (#, п) + k2и * (#,7) = -8(£, ф, Г/Є Т.

и

к ,1

и

N

к ,1 =0

Г

и

о

к .1

Фундаментальное решение для уравнения Гельмгольца в двумерном однородном пространстве известно и равно

и*

где г =

>/[х(п) - х(#)]2 + [у(п) - у(#)]

Н,(1) (Кг) = ./о(Кг) + /Т0(Кг) - функция Ханкеля первого рода и нулевого порядка, где J0 - функция Бесселя первого рода нулевого порядка, У0 - функция Неймана нулевого порядка.

Для построения граничного интегрального уравнения для рассеянного поля и его нормальной производной в зонах р и Т3 воспользуемся теоремой Грина в следующем виде:

| ^2и* (#, г;) + к V (#, 7)] и(п)ё ¥

¥

= -[ КпК (#,п)ёг +

и

где у(п) =

ди (#,п) дп'

ди(п) дп'

ё г +1 /р (7)и* (#,п)ё ¥,

нормальные производные поля,

под Г' принимаются бесконечные границы у = а и у = 0 для зон р и Т3 соответственно.

Подставляя равенство (16) в уравнение (18) и осуществляя предельный переход точки наблюдения % из внутренней в граничную, найдем

с(#)и(#) = -| /р (п)и* (#,п)а т +

т

д .

+| у(ц)и* (#,п)а г -1 и (п) и да,л а г.

Это уравнение обеспечивает связь между полем и и его нормальной производной V на границе Г'. Функция с(#) в уравнении (19) равна:

с

2п

1 ™ I

- иг, (х, ) =2 h К (хт ) 2

'к0а0

где ф - внутренней угол кусочно-линейной границы в точке #. Первое слагаемое в правой части (19) является полем, создаваемым источником в свободном пространстве, и может быть обозначено как падающее поле ир.

Таким образом, уравнение (19) записывается в следующем виде:

с(#)и(#) = | Кп)и* (#, п)ёг -

-Г и(п) ди (#,П) ё г + ир (#).

•* Яп

дп

Поле и и его производные V в случае дифракционной решетки являются квазипериодическими функциями [11, 13-15]:

и(X, у) = и(х, у)е' К X, у) = я X, у)е'"

где и( х, у) и V( х, у) - периодические по х функции с периодом а.

Интегрирование по бесконечной границе Г' в (21) можно заменить интегрированием по границам Г и Г3:

2и( х0> у) =

Г 2 Vх, у)е*0а(х+пё)и (х + пё - х0) ёх -

г п=-да да

-[ 2 и(X,у)е"

(х+пё) ди (х + пё - х0) ёх +

дп'

+и'п (х0, у), х є [0, ё),

где Г = Гь Г3, у = а и у = 0 для границ Г1 и Г3 соответственно.

Подставляя вместо функции комплексной амплитуды в уравнении (23) на границе ее аппроксимацию базисными кусочно-линейными функциями (8) и (9), получим соотношение (24).

[ аТ' (хт + ^~)иТ (хт + пё - х, + кц) ёц

-иг (хт) 2

"к0а0

Г г/ , ЧдиТ'( хт +пё - X + ^),

[  ' (хт + ^------------------—------------ёП

дп'

>+(х,),

2

і

г

Г

г'

г' п=-да

т=1

п=-да

-1

п=-да

-1

где 5 = [1, ^-1], I = 1, 3.

Граничные условия на поле и его производные на границах Г2 и Г4 являются периодическими вида:

иГ = иГ е*0^, vГ = VГ е^.

Г 4 Г 2 Г 4 Г 2

Разложения (8) и (9) и граничные условия (25) позволяют записать для коэффициентов поля и его производных на границах Г2 и Г4:

и

где 5 = [1, N1. Элементы матриц D и G, соответствующие границам Г2 и Г4, записываются как:

а25 = 1, а45 = а2У*0“0а, Я25 = 1, Я4 = g2vw, е25 = 0, е4 = 0,5 = [1, Ny\.

Соотношения (24) могут быть представлены в матричном виде (30) с элементами матриц D, G и Е вида (27), (28) и (29). Бесконечные ряды в (28) и (29) аппроксимируются конечными суммами с относительной погрешностью менее 10-5.

dr = -h

*ОаО

gi

= h

'кОаО

І du* (Xm + nd - Xs + hn),

j аг- (Xm + hn)--------^--------—----------------dn

dn'

eS = -s,

i 2 s’

| or‘ (xm + hn)u* (xm + nd - xs + hn) dn

m, s = [І, Nr-І], i = І, 3. [D]ur + [G]vr = [E] Un.

Интегралы в равенствах (28) и (29) могут быть оценены численно. Объединяя системы уравнений (15) и (30), получим замкнутую систему линейных алгебраических уравнений для решения задачи дифракции плоской волны на периодическом двумерном микрообъекте (31).

[aq,q] [at,q] О] ua

[Ао,г] [Аг,г] [B] «г

. О [D] [G]_ . Уг_

"[CQ,Q] [Сг,а] О]  fn‘

[СП,г] ^[Сг,г] О fг

О О [e]_ _<.

где подматрица An,n размерностью (N - M)*(N - M) включает в себя коэффициенты соотношений поля во внутренних узлах сети разбиения, подматрицы АПГ и АГ,П размерностью (N - M) *M и M*(N - M) соответственно включают в себя коэффициенты связи поля в граничных и внутренних узлах, подматрица АГ, Г размером M*M включает в себя коэффициенты связи поля в граничных узлах, подматрица В размером M*M включает в себя коэффициенты соотношений между производными поля на границе и полем во внутренних узлах сети разбиения, подматрица D размером M*M включает в себя коэффициенты связи поля свободного пространства в граничных узлах, подматрица G размером M*M включает в себя коэффициенты соотношений между производными поля на границе и полем свободного пространства во внутренних узлах сети разбиения. Подматрица Cn,n размерностью (N - M)*(N - M) включает в себя коэффициенты соотношений источников поля во внутренних узлах сети разбиения, подматрицы СП,Г и СГ,П размерностью (N - M) *M и M*(N - M) соответственно включают в себя коэффициенты соотношений источников поля в граничных и внутренних узлах, подматрица СГ, Г размером M*M включает в себя коэффициенты соотношений источников поля в граничных узлах. Подматрица E имеет размер M*M Вектора un и иГ - вектора напряженности поля во внутренних и граничных узлах сети, vr - вектора нормальных производных поля в граничных узлах сети. Вектора fn и fr - вектора источников поля во внутренних и граничных узлах сети, и” - вектор напряженности внешнего падающего поля в граничных узлах сети. Таким образом, система (31) имеет N + M уравнений и столько же неизвестных.

Поле в областях ^ и Т3 может быть представлено в виде разложений Рэлея (суперпозиций плоских волн). В зоне г-компоненты полей имеют вид:

и (х, у) = ехр (/^ (а х - ву)) +

”

+ 2 Кп ЄХР (г10 (апх + РпУ)) ,

где ап =V^ sin

Функция и (x, у) представляет компоненту Ez (x, у) электрического поля для случая ТЕ-поляризации и компоненту Hz (x, у) магнитного поля для случая ТМ-поляризации.

В зоне Т3 z-компоненты имеют вид:

и

п=-да

где вп = \М -ап2 .

Разложения Рэлея (32), (34) являются решениями уравнения Г ельмгольца

Дu + к u = О ,

при k2 = k02e1 и k2 = к^є3, соответственно.

После определения значений поля в области О и его производных на границах Г1 и Г3 нормированные интенсивности отраженных и прошедших порядков рассчитываются следующим образом [1, 5, 6]:

IR = \R"\

в

в

IT =,

2 Рп_

во

I IR + I IT = І

tR In I2 вt tT і''і І7- і2 вn

It = К| в її W в

(

I IR + I ITT = І

для ТЕ- и ТМ-поляризованных волн соответственно. Выражения (37), (38) определяют, какая часть энергии падающей волны перейдет к п-му порядку дифракции. Заметим, что интенсивности порядков распространения существуют для непоглощающих материалов, т.е. для 1т.§£1 = 0. Множества и и и2

в (37) и (38) являются множествами индексов, соответствующих отраженным и прошедшим порядкам дифракции:

"=-да

"=-да

"=-да

2

V'^l

' є Z: — < 11-, если Іт^Є = 0

Є

0,

если Іт^Іе > 0

п є Z

0,

Для определения коэффициентов отражения Яп и пропускания Гя воспользуемся дискретным преобразованием Фурье:

1 N-1

Кп = 2[и

N к=0

1 N-1

Тп = Т?Е и

N к=0

Эти коэффициенты описывают амплитуду и фазовый сдвиг распространяющихся плоских волн. Более точно, модули |Лп| и |Тп| - это амплитуды п-го отраженного и прошедшего порядков и а^ [Яп/\Яп |] и агё [Тп/\Тп|] - фазовый сдвиг. Коэффициенты с п г Ц п й Ц3 описывают затухающие волны.

Результаты численного моделирования

Рассмотрим в качестве примера дифракцию плоской волны с длиной волны Л0 = 0,6 мкм на бинарной диэлектрической дифракционной решетке с коэффициентом заполнения 0,25 с толщиной 0,24 мкм. Период решеток менялся от 0,2 мкм до 2 мкм. Соответственно при моделировании ПМКЭ-ГЭ длина сегмента сети покрытия к менялась от Л0 / 300 до Л0 / 30 .

Рассмотрим пример. Плоская волна падает по нормали из воздуха (е1 = 1) на решетку (е3 = 2,25). На рис. 1 и 2 представлены зависимости эффективностей I? и

/^ порядков дифракции ТЕ- и ТМ-поляризованной волны соответственно от периода решетки.

Зависимости отклонений интенсивностей А/? и

д/Т0 (а/ = \/?

нулевых порядков ди-

фракции, рассчитанных RCWA и ПМКЭ-ГЭ методами, ТЕ- и ТМ-поляризованных волн от длины сегмента сети покрытия h метода ПМКЭ-ГЭ приведены на рис. 3. Трудно отметить явную зависимость отклонения для h < Л0 /30 . Равномерная норма отклонения интенсивностей составляет 4-10-3 и 2-10-3 для ТЕ- и ТМ-поляризаций соответственно.

В следующем примере плоская волна падает по нормали из подложки решетки (£1 = 2,25) в воздух £ = 1). На рис. 4 и 5 представлены зависимости интенсивностей и II порядков дифракции ТЕ- и ТМ-поляризованной волны от периода решетки.

Зависимости отклонений интенсивностей Д1^ и Д 1Т0 порядков дифракции ТЕ- и ТМ-поляризованных волн от длины сегмента сети покрытия h приведены на рис. 6. Равномерная норма отклонения интенсивностей составляет 7-10-3 и 8-10-3 для ТЕ- и ТМ-поляризаций соответственно.

— < 11-, если Іт^Є = 0

є 1

если Іт у[е > 0

а)

Период, мкм

б)

Период, мкм

Рис. 2. Распределение эффективности нулевых порядков дифракции ТЕ-поляризованной волны: а) I^ , б) 1Т

а)

Период, мкм

б)

Период, мкм

Рис. 3. Распределение эффективности нулевых порядков дифракции ТМ-поляризованной волны: а) I?, б)

)

а)

ДR0

ДТ0

Длина сегмента сети, мкм

дгаэ

дто

б)

Длина сегмента сети, мкм

Рис. 4. Зависимость отклонений эффективности нулевых порядков от параметра к: а) ТЕ-волны; б) ТМ-волны

а)

Период, мкм

б)

Период, мкм

Рис. 5. Распределение эффективности нулевых порядков дифракции ТЕ-поляризованной волны:

а) К,б) 1Т

а)

Период, мкм

б)

Период, мкм

Рис. 6. Распределение эффективности нулевых порядков дифракции ТМ-поляризованной волны:

а) 1о • б) 1т

ДR0

ДТ0

а)

Длина сегмента сети, мкм

ДR0

ДТ0

б)

Длина сегмента сети, мкм

Рис. 7. Зависимость отклонений эффективности нулевых порядков от параметра к: а) ТЕ-волны; б) ТМ-волны

Таким образом, сравнение результатов, полученных разработанным ПМКЭ-ГЭ методом и RCWA методом, показало их хорошее соответствие.

Заключение

В данной работе представлен объединенный метод конечных элементов Галеркина и граничных элементов для решения задач дифракции плоской электромагнитной волны на периодических неоднородных двумерных (цилиндрических) объектах микрооптики, период которых сравним с длиной волны.

Проведено исследование зависимости относительной погрешности результатов моделирования от отношения длины волны падающего света к длине сегмента сети покрытия.

Благодарн ости

Работа выполнена при поддержке российско-американской программы «Фундаментальные исследования и высшее образование» (грант CRDF RUX0-014-SA-06), гранта РФФИ 08-07-99007-р_офи, Фонда содействия отечественной науке и гранта Президента РФ поддержки ведущих научных школ (НШ-3086.2008.9).

Литература

1. Electromagnetic Theory of Gratings: Topics in current physics / Ed. by R.Petit // N.Y.: Springer-Verlag. - 1980.

2. Moharam, M.G. Rigorous coupled wave analysis of planar grating diffraction / M.G. Moharam [and others] // J. Opt. Soc. Amer. 71. -1981. - P. 811-818.

3. Chandezon, J. A new theoretical method for diffraction gratings and its applications / J. Chandezon [and others] // J. Opt. (Paris), 11. - 1980. - P. 235-241.

4. Urbach, H.P. Convergence of the Galerkin method for two-dimensional electromagnetic problems / H.P. Urbach // SIAM J. Numer. Anal., 28. - 1991. - P. 697-710.

5. Bao, G. Finite element approximation of time harmonic waves in periodic structures / G. Bao [and others] // SIAM J. Numer. Anal., 32. - 1995. - P. 1155-1169.

6. Elschner, J. Finite element solution of conical diffraction problems / J. Elschner [and others] // Advances in Comp. Math., 16. - 2002. - P. 139-156.

7. Prather, D.W. Boundary integral methods to the analysis of diffractive optical elements / D.W. Prather [and others] // J. Opt. Soc. Am. A, 14. - 1997. - P. 34-42.

8. Yee, K. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media / K. Yee // IEEE Trans. Antennas Propag. AP-14. -1966. - P. 302-307.

9. Saj, W.M. FDTD simulations of 2D plasmon waveguide on silver nanorods in hexagonal lattice / W. M. Saj // Opt. Express 13. - 2005. - P. 4818-4827.

10. GSolver, rigorous diffraction grating analysis [Электронный ресурс]. http://www.gsolver.com. - Grating Solver Development Company.

11. Методы компьютерной оптики / под ред. В.А. Сой-фера. - М.: Физматлит, 2000. - 688 с.

12. Неганов, В.А. Линейная макроскопическая электродинамика / В.А. Неганов, С.Б. Раевский, Г.П. Яровой. -М.: Радио и Связь, 2000. - 509 с.

13. Methods For Computer Design of Diffractive Optical Elements / Ed. by V.A. Soifer // N.Y.: Wiley-Inter-science Publication John Wiley & Sons, Inc. - 2002.

14. Moharam, M.G. Formulation for Stable and Efficient Implementation of the Rigorous Coupled-Wave Analysis of Binary Gratings / M.G. Moharam [and others] // J. Opt. Soc. Amer. 12 (5). - 1995. - P. 1068-1076.

15. Нестеренко, Д.В. Анализ дифракции света на элементах цилиндрической микрооптики объединенным методом конечных элементов Галеркина и граничных элементов / Д.В. Нестеренко, В.В. Котляр // Компьютерная оптика, 32. - N1. - 2008. - C. 23-28.