CC BY
14
Подходы к алгоритмизации строгого метода связанных волн
А.И. Антонов', Л.А. Васин', Г.И. Грейсух1 'Пензенский государственный университет архитектуры и строительства, Пенза, Россия
Аннотация
Выведены коэффициенты Фурье в разложении в ряд диэлектрической проницаемости линейно-пилообразного рельефа с положительным тангенсом угла наклона рабочей поверхности и с вертикальным обратным скатом, необходимые для алгоритмизации и реализации строгого метода связанных волн. Описаны и сравнены по устойчивости и эффективности подходы к реализации строгого метода связанных волн: подход целенаправленного преобразования матрицы пропускания и подход Гауссовых сокращений. Сделан вывод, что оптимальным подходом для определения дифракционной эффективности в случае пилообразной микроструктуры, глубина рельефа которой одного порядка с пространственным периодом, является подход целенаправленного преобразования матрицы пропускания.
Ключевые слова: дифракционная эффективность, диэлектрическая проницаемость, уравнения Максвелла, строгий метод связанных волн.
Цитирование: Антонов, А.И. Подходы к алгоритмизации строгого метода связанных волн / А.И. Антонов, Л.А. Васин, Г.И. Грейсух // Компьютерная оптика. - 2019. - Т. 43, № 2. -С. 209-219. - БО!: 10.18287/2412-6179-2019-43-2-209-219.
Введение
Согласно представлениям современной физики, электромагнитное излучение исчерпывающе описывается системой уравнений Максвелла. Однако решение этой системы для определения дифракционной эффективности (ДЭ) различных дифракционных оптических элементов (ДОЭ) без дополнительных упрощений представляет собой, как правило, весьма сложную задачу.
Скалярная теория дифракции (СТД) позволяет рассчитывать ДЭ с меньшими затратами времени и труда по сравнению с электромагнитной теорией дифракции (ЭТД). При этом в ряде случаев наблюдается хорошее совпадение результатов, даваемых СТД и ЭТД. Однако для ДОЭ с относительно низким (меньше десяти) значением отношения периода дифракционной микроструктуры к высоте рельефа результаты, получаемые в рамках СТД, значительно превышают значения, даваемые ЭТД [1, 2, 3, 4]. В таких случаях достоверность результатов может быть гарантирована только путем решения уравнений Максвелла.
Существует несколько методов решения указанных уравнений, и в каждом из них используются определенные допущения. Это относится как к аналитическим методам (например, к строгому методу связанных волн (RCWA)), так и к численным методам решения дифференциальных уравнений (например, к методу конечных разностей во временной области (ГОТО)).
Алгоритмы реализации каждого метода также могут различаться. В данной статье описаны два основных подхода к реализации RCWA-метода. Исторически первым был подход, предполагавший прямое решение системы уравнений, описывающих граничные условия для уравнений Максвелла [5]. Данный подход допускает сокращение количества неизвестных методом Гауссовых преобразований. В результате в системе уравнений, описывающих граничные условия, останутся только амплитуды отраженных и прошедших порядков дифракции. Указанное сокращение можно рассматривать как определенную модерниза-
цию данного подхода, ориентированную на компьютеры с невысокой производительностью - подход Гауссовых сокращений (ГС-подход).
Второй подход основан на целенаправленном преобразовании матрицы пропускания (ЦПМП). Он предполагает решение системы дифференциальных уравнений второго порядка, приводящей к системе линейных алгебраических уравнений, описывающих граничные условия, с последующим преобразованием и упрощением этой системы, исключающей неустойчивость решения благодаря хорошей обусловленности обращенных матриц [6, 7]. Количество неизвестных констант в системе уравнений граничных условий в данном подходе больше, чем в соответствующей системе граничных условий первого подхода за счет решения системы дифференциальных уравнений второго порядка. Однако все неизвестные константы исключаются в процессе получения итоговой матрицы пропускания.
Общими положениями RCWA-метода являются следующие допущения. Периодическая дифракционная микроструктура представляется совокупностью последовательно расположенных вдоль одной из трех взаимноперпендикулярных осей диэлектрических слоев определенной толщины (например, вдоль оси 2, как это показано на рис. 1). При этом диэлектрическая проницаемость (ДП) каждого слоя одинакова по всей толщине йп и периодична вдоль двух других осей еп (х,у) = еп (х +Л х,у +Лу), где п - номер слоя, Л х, Л у - периоды вдоль соответствующих осей. В результате ДП каждого слоя можно представить в виде двумерного ряда Фурье, при этом точность моделирования рельефа поверхностной микроструктуры определяется количеством слоев, на которые разбивается рельеф.
Далее составляющие векторов электрической и магнитной напряженностей внутри слоев и в двух средах, из которой излучение падает на микроструктуру и в которую входит после микроструктуры (области ! и !! на рис. 1), также представляются в виде
одномерных либо двумерных рядов Фурье. В области I распространяются отраженные порядки дифракции, а в области II - прошедшие. С использованием выражений для электрических и магнитных напряженно -стей внутри каждого слоя и условия непрерывности составляется система уравнений граничных условий, которая решается относительно амплитуд отраженных или прошедших порядков дифракции. Здесь заметим, что в настоящей статье с целью уменьшения объема представляемых математических выкладок, но без потери общности, рассматриваются наиболее часто используемые киноформные микроструктуры, которые являются периодичными только вдоль одной оси и у которых пространственный период одного порядка с глубиной рельефа. Что касается количества слоев, на которые разбивается рельеф Ы, то оно выбирается из условия, обеспечивающего практически предельную для анализируемых подходов точность вычисления ДЭ при минимальных затратах времени на расчет: М»ШИ / X. Такое число слоев ниже будем называть оптимальным.
О
Область I
п=1 ' ' 1 1 , —— Г—
п=т - ц г1-ц — 11 11 ' А*
. 1 1
и=ЛГ 1 1
1 1
гБ
к
Область II
Рис. 1. Пример разбиения одного периода поверхностного рельефа с глубиной И микроструктуры ДОЭ на N слоёв: кпс - волновой вектор падающей на микроструктуру плоской волны
Завершая введение, отметим, что, наряду с вышеперечисленными подходами, существует и так называемый подход матриц рассеяния [8]. Данный подход предполагает многочисленные матричные операции, в том числе такую, как матричное произведение Рэд-хеффера [9]. В результате при расчете ДЭ микроструктур с периодичностью только вдоль одной оси и с относительно небольшим, но оптимальным количеством слоев RCWA-метод, базирующийся на подходе матриц рассеяния, значительно уступает в скорости методу, базирующемуся на ЦПМП-подходе.
В данной статье рассматриваются основные положения ГС-подхода и подхода ЦПМП. В оригинальных статьях [5 - 7] не описаны приемы численной реализации данных подходов, особенно это относится к работе [5], использование результатов которой невозможно без преобразования особым способом матрицы граничных условий. В работе же [7], например, не представлен алгоритм, по которому можно получить итоговую матрицу прошедших дифракционных порядков, являющуюся множественным матричным произведением (см. параграф 4). Также в [5 - 7] не описаны способы получения коэффициентов Фурье в разложении ДП в ряд для любой микроструктуры и, в частности, не получены выражения коэффициентов
для пилообразной микроструктуры, анализируемой в этих работах и в настоящей статье.
1. Определение коэффициентов Фурье в разложении ДП
В рассматриваемом случае периодичности вдоль оси Х ДП п-го слоя имеет вид еп (х) = еп (х+Л) и может быть описана рядом Фурье:
-(х )=!
0]аКх
1 Л
/е-п (х к1акх
(1)
(2)
где а - номер гармоники, ] - мнимая единица, /Е п (х) -функция, задающая значение ДП в каждой точке слоя, К = 2п / Л.
Функция /гп (х) является бинарной и принимает в пределах периода Л только два значения - е1 и е2. В зависимости от формы рельефа интеграл в (2) может быть представлен в виде суммы интегралов по однородным участкам периода, и каждый из таких интегралов будет содержать множитель е1 или е2.
Пусть в п-м слое рельефа микроструктуры участок с ДП, равной е2, находится внутри периода между точками с координатами /1п и /2 п , как показано на рис. 2.
О /;,„ 12,п А X
21 £] е2 е1
Рис. 2. Распределение ДП вдоль оси X В этом случае подынтегральную функцию /¡¡,п(х) можно записать как
А,п (х) =
е,
х 6 Е1,п
х 6 Е2,п '
где Е1п - участки в пределах периода Л с ДП, равной е,, а Е2п - участок в пределах периода с ДП, равной е2. Тогда после соответствующих преобразований интеграла (2) получим:
£ а,п
(е2 -е,)— | е-'аКхах аф 0
Л
(3)
е2Л + е, (1 - /я) а = 0,
где /п = (12,п - /,,„) / Л - доля периода Л с ДП, равной е2, в слое с номером п. Выражение (3) используется в качестве коэффициента Фурье в подходе ЦПМП. Для ГС-подхода [5] ряд (1) преобразуется в
еп (х) = е, +(е2 - е, ^е^п • е]аК, (4)
где
е =
°ап
, Л
Лп (х)• е-^йх,
л
а.п
е
'2,п
а.
Л (х) =
О х е Е^ х еЕ9 п
Очевидно, что для слоя с двумя значениями ДП
р = -
'2,п
| е-]оКх¿х.
Все возможные распределения двух ДП в слое будут различаться коэффициентами £ап, зависящими от границ интегрирования в (3), при этом сами эти границы для каждого слоя определяются формой рельефа микроструктуры.
Наиболее часто используемая киноформная микроструктура имеет, как правило, линейно-пилообразный рельеф с положительным тангенсом угла наклона рабочей поверхности и с вертикальным (параллельным оси 2) обратным скатом (рис. 3). Пусть число слоев, на которые разбивается рельеф, равно N. Правая вертикальная граница /2,и завершает период Л, а левая 11п является функцией номера слоя. Легко показать, что
кп = > - п).
12.т=Л X
п = 1 —3 г 1 ? С?
п ! Ь'
у ✓ у ! /
п — т Ъ''
и =N-1 Е/ у у ✓ ✓ ✓ ✓ е 2
п=Ы V
✓
А г
А
N
Рис. 3. Пример изменения от слоя к слою ширины участков с ДП, равной е2; к - глубина рельефа
Тогда для рассматриваемого пилообразного рельефа коэффициенты Фурье в выражении (1) будут иметь вид
р =
( - £2 )
Ле'
- ] 2па
(еЛ2 ™}п -1)
а Ф О
2па
I + £1 (1 -1) а = О,
(5)
гдеЛп = п / N - доля р2 в слое п. Отметим, что £0,п является средним значением ДП слоя.
2. ЦПМП-подход для случая ТЕ-поляризации
Применив к уравнениям Максвелла, описывающим электромагнитное поле в изотропной среде, метод комплексных амплитуд, получим систему уравнений
'дИ2 дИу . дЕу Ж2 .
ду дг дг ду
дИх дИг . дЕг дЕх . „ ...
—---— = ¡о^еЕ, —---— = ,/юЦоЦИу, (6)
дг дх дх дг
дИу дИх , „ дЕх дЕу .
—---— = ;юеорЕг , —---— = ¡^0^ ,
дх ду ду дх
где Их, Иу, Иг, Ех, Еу, Ег - составляющие векторов на-пряженностей магнитного и электрического полей; д, £ - значения магнитной и диэлектрической прони-цаемостей рассматриваемой среды.
В RCWA-методе составляющие векторов электрической и магнитной напряженностей внутри слоя представляются в виде рядов Фурье.
В рассматриваемом случае периодичности только вдоль одной оси указанные ряды можно записать как [10]
Ех = £ Бх, (г)е
да
Еу = £ Бу; (г)е
/=-да да
Ег = £ Б, (г)е
- ¡К,
- ¡кх,,
Их = Л - £ их,, (гу*,
\ Д0 ,=-да
(7)
- ]кх„
Иу = Ц — £ иу, (г) д0
Иг = Л£ иг,, (г)е
V Д 0 1=-да
х-составляющая волнового вектора волны, дифрагированной в 1-й порядок, (кх/) находится как следствие теоремы Флоке [9] и периодичности ДП вдоль оси X:
кх, = кхпс -/Кх, / е Z; (8)
где кхМс - х-составляющая волнового вектора излучения, входящего в среду; Ъ - множество целых чисел. При реальных расчетах ряды Фурье (7), а следовательно, и количество учитываемых дифракционных порядков, естественно, ограничиваются: /тт < / < /тах.
Пусть падающая на пилообразную микроструктуру плоская электромагнитная волна имеет ТЕ-поляризацию, а магнитная проницаемость всех слоев д = 1. Тогда, подставляя (7) в (6), учитывая, что ю = к0 / (£0д0)1/2, и объединяя суммирования в (1) и (7), получаем матричные уравнения [10]:
д
и п,х + ¡Кх и п, г = Е п8 п
д(г)
__8 = и
д(2') пу Ип
Кх8п,у = ¡ип,г ,
(9)
где индекс п показывает принадлежность матриц к п-му слою (см. рис. 1), Кх - диагональная матрица с элементами к'х,, = кх,,/ к0, х'= к0г, к0 = 2п / X (X - длина волны падающего на микроструктуру излучения); 8„у, И пх, И п,г - матрицы-столбцы элементов Бпу/ (г'), П„хХ,г (г') и Пп,г,,(г') с индексами / от /'тт до /'тах; Е -двумерная матрица с элементами £/-и (/ - номер строки, и - номер столбца); М - матрица магнитных про-
! >
ницаемостей, аналогичная по форме матрице E. Для случая TM-поляризации можно получить систему матричных уравнений, аналогичную (9).
Решив систему уравнений (9) относительно матрицы-столбца Sny, получим выражения для амплитуд Sn,y,i (z ), зависящих от неизвестных констант и элементов матрицы собственных векторов WE,n матрицы BE,n = Kx2- En [6, 12]. В RCWA-методе напряженности электрических и магнитных полей в областях I и II также представляются в виде рядов Фурье; в случае TE-поляризации в нуль не обращаются только y-составляющие векторов напряженности электрического поля Ey в областях I и II. Коэффициенты Фурье для области I (R) представляют собой амплитуды отраженных порядков дифракции, для области II (T) -амплитуды прошедших порядков дифракции [6]. х-составляющие векторов напряженности магнитного поля Hx можно определить из выражений для Ey и системы уравнений (6). Амплитуды х-составляющих векторов магнитной напряженности внутри каждого слоя Unyi(z') находятся из уравнений системы (9).
Компоненты вектора электрической напряженности, касательные к границе раздела двух сред, всегда непрерывны, и, поскольку рассматривается случай без поверхностных токов, на данной границе также непрерывны компоненты вектора магнитной напряженности. Тогда, последовательно приравнивая Ey каждой пары граничащих друг с другом слоев, включая области I и II, и проделывая аналогичную операцию для Hx, можно получить полную систему граничных условий, содержащую неизвестные константы и амплитуды Ri и Ti. Можно произвести исключение неизвестных констант. Однако операция непосредственного решения получившейся системы относительно Ri и Ti требует обращения матриц, содержащих экспоненциальные элементы exp(- qn^mk0dn) (здесь qnm - квадратный корень m-го собственного числа матрицы BEn, dn - толщина n-го слоя (см. рис. 1)). Количество операций обращения будет равно количеству слоев. Ввиду своей малости экспоненциальные элементы могут сделать эту матрицу почти сингулярной. В результате даже при аналитическом нахождении обращенной матрицы сингулярным разложением [15], из-за усечений с конечной точностью в числовом представлении получившихся больших чисел, найденные матрицы будут ошибочными, что приведет к неверному результату и неустойчивости численного решения всей системы.
Можно сделать такие преобразования сокращенной системы уравнений граничных условий, чтобы избежать обращения матриц, содержащих экспоненциальные элементы exp(- q„mk0dn) [7]. Тогда окончательно получим удобное для вычислений матричное уравнение:
5,о
n j cos (0 )5i-0
-I f1" " R"
_ jke ,i g1 _ _T1 _
(1о)
где Я - матрица-столбец элементов Я,-; Т - матрица-столбец элементов Т; 5,0 - матрица-столбец элемен-
тов дельт Кронекера [14] с целыми /тш < /' < /тах; I -единичная матрица размерности (5 - учитываемое количество дифракционных порядков); - диагональная матрица элементов к^ / к0, q = 1,11;
kz,q,i =>/(kq )2 - (kx,i )2 , q
= I, II;
q - номер области; (к^ = (2^ / X) и (кп) = (2лпп / X) - квадраты модулей волновых векторов отраженных (область I) и прошедших (область II) дифрагированных волн соответственно; ^ - показатель преломления среды в области q;
Wen (I + X„b„a-'Xn)
gn _ VE,n ( — Xnbиа- ■X,)_
a„ " We ,n WE ,n ' -1 fn+1
b „ _ _ Ve n -Ve n _ g n+1 _
причем ^+1 = I, §N+1 = /Ке,п, Уб,п = WE,ИQИ, - диагональная матрица элементов qn,m, Хп - диагональная матрица элементов ехр(- q„,mkodn), Тп = ап^Хп-иТ- [7]. Следовательно,
Т = а-/ХN ■■■а
- Хп ■ ■ .Т1 . (11)
Таким образом, решая систему линейных алгебраических уравнений (10), можно найти все Я, а подставляя матрицу-столбец Т1 из решения (10) в выражение (11), - все Т,-.
3. ГС-подход для случая ТЕ-поляризации
Рассмотрим случай ТЕ-поляризации с ранее выбранными характеристиками областей I и II Из системы (6) получим уравнение:
d En y d En
+ k0 8nEn, y =
(12)
дх2 дг2
Представим касательную составляющую вектора электрической напряженности в виде ряда:
да , ~ ,
- ]\кх ,,х+кг ,пг)
En,y ^ ' y(z)e
(13)
В отличие от (7) в выражении (13) выделена экспонента ехр(-/кгпг), зависящая от г и содержащая г-составляющую волнового вектора внутри п-го слоя кгп [5]. Поскольку ДП не периодична в направлении г, волновые векторы всех дифракционных порядков имеют одну и ту же компоненту волнового вектора в этом направлении, то есть г-составляющую нулевого (не измененного по направлению) дифракционного поряди
Подставив ряды (4) и (13) в уравнение (12) и сделав замену второй производной от амплитуд 8п,у,(г), можно получить систему дифференциальных уравнений первого порядка, с помощью которой определяются сами амплитуды £пу,,(г) [5] Произведя операции с граничными условиями, аналогичные описанным в ЦПМП-подходе, получим систему уравнений граничных условий ГС-подхода. Пусть
0 I
ВЕп = . _. , (14)
i
где
A n —
'imin an,-l • • • an,-(s-l)
an,l ' ' an,-l :
: an,i '. an,-i
an,s-l ''' an,1 bimax
an
an,a — -ко (nii - n2 )e^; bt — -iK(2kxjK - iK);
— 2 jkz,n ; kn — \ ^ O n
^ ^ 0,n ^^ ; kz,n — ^kn kX
и W'E,„ - матрица собственных векторов матрицы Б'Е,„, ограниченная первыми 5 строками с (р, т)-элементами w„p,m (1 <р < 5, 1 < т < 25); р пт - т-е собственное число матрицы Б'Е,п.
Тогда полученная система граничных условий запишется в матричной форме:
H
Cf R T
— v,
(15)
где
н — [ч % z]; * — [Ф1
ф
n
ф Л
ф
n
; % —
" -I" " Z 4 " §,о
-KI ; Z — -I ; v — -KI 5,о
_ Z _ K ii _ _ Z5 _
ф1, п ф2,п фз,п ф4,п Z2,n
ср - матрица-столбец неизвестных констант с количеством элементов 25^ Zl,„ - нулевая матрица размером 25(„+1)*25; Z2,„ - нулевая матрица размером 25^+1)*25; Z3 - нулевая матрица размером 25^5; Z4 - нулевая матрица размером 25^5; Z5 - нулевая матрица-столбец размером из 25N элементов; К': -диагональная матрица элементов ¡к^,; К'П - диагональная матрица элементов ¡к^ц; ф1,„, ф2,„, фз,„, ф4,„ -(р,т)-матрицы размером 5*25 с элементами соответственно
H5ln+1 Wn, PJ( - jkzn )
(-l)8l n +1 Wn,pm (pm,n - jkz,n ) e
Pm.n - jkzn )Dn-
(множитель (-1)81,„ обеспечивает снятие знака «минус» при п = 1),
( - jkz n )
Pmn - jk z ,n )Dn ,
Wn,p,m (m,n - jkz,n )
n N
Dn — £ dr, Dn — £ dr — h
(см. рис. 1). Матрица Н имеет размер 25^+1)*25^+1). В правой части системы (15) взята
матрица -К'До, поскольку только один ее элемент не обращается в нуль (при / = 0, а кгП,0 = кг/„с).
Система (15) может решаться относительно Я и Т непосредственно либо с предварительным сокращением всех элементов матрицы ср (ГС-подход). Методика сокращения системы (15) преобразованиями Гаусса будет описана в параграфе 4.
4. Алгоритмы расчета ДЭ ЯСЖ4-методом при использовании ЦПМП- и ГС-подходов
Выражение для ДЭ в общем случае можно получить, используя сонаправленность волнового вектора и вектора Пойтинга электромагнитной волны. В результате для ДЭ отраженных и прошедших порядков дифракции в случае одномерных ДОЭ будет получено [8]:
DEr , — RR* Re
(К, Л
DET f — TT* Re
k
у ^z,inc у (k Л
k
у ^z,inc у
(16)
Для анализа и сравнения двух подходов к определению амплитуд, прошедших и отраженных дифракционных порядков, а именно подходов ЦПМП и ГС, была рассмотрена одномерная киноформная микроструктура с линейно-пилообразным рельефом, с положительным тангенсом угла наклона рабочей поверхности и с вертикальным (параллельным оси Z) обратным скатом (рис. 3). Коэффициенты Фурье описывались выражениями (7). В качестве исходных параметров системы были взяты: показатели преломления областей I и II (nI и nII), длина падающей на рельеф волны X, угол ее падения (предполагалось, что плоскость падения совпадает с плоскостью X- Z) 9, глубина рельефа h, отношение периода решетки к глубине рельефа (P = Л / h), число учитываемых дифракционных порядков s, число слоев, на которые разбивается рельеф N. Был выбран случай TE-поляризации.
Первоначально в таких работах, как [1 - 4], для определения ДЭ RCWA-методом использовался закрытый алгоритм, разработанный сотрудниками Института общей физики им. А.М. Прохорова Российской академии наук (ИОФ РАН) [15]. Также существуют коммерческие продукты, рассчитывающие ДЭ RCWA-методом (например GD-Calc [16]) и использующие для реализации среду MATLAB [17].
Рассмотрим предлагаемые авторами статьи численные реализации ЦПМП- и ГС-подходов. Для написания алгоритмов выбрана среда Wolfram Mathematica 11.2 [18]. Далее описаны общие положения реализации этих алгоритмов без обозначения внутреннего синтаксиса программы.
Период микроструктуры во всех расчетах определяется через введенный выше параметр P, а каждый слой микроструктуры считается состоящим из двух участков с ДП £i и ец, т.е. теми же, что и в областях I и II.
С помощью внутренних инструментов системы Wolfram Mathematica возможно создавать многомерные матрицы с любой очередностью индексов или аргументов каждого элемента-функции этой матрицы. Таким образом, создаются не только непосредственно двумерные матрицы En и BE,n, которые, в свою очередь, являются функциями от номера слоя n, а следовательно, различны и считаются отдельно для каждого слоя, но и списки целых матриц, например, В - список матриц BE,n всех слоев, что вводится для оптимизации вычислений. Напомним, что En - матрица элементов £,_ди (eI = nI2, eI = nII2 ), где i - номер строки, p - номер столбца элемента матрицы; BE,n = Kx2- En. Здесь Kx2 - диагональная матрица элементов
Ж =
v kt) =(nisin (e)-' Л).
Аналогично B строятся матрицы W, Q и Х, как, соответственно, список матриц WE,n собственных векторов BE,n для каждого слоя, список диагональных матриц Qn с элементами - квадратными корнями собственных чисел BE,n, список диагональных матриц Xn c элементами exp(-2nqnmd / (XN)), где qn,m - m-е собственное число BEn. Далее строится список V скалярных произведений WE,n и Qn.
Особенности построения списков в среде Wolfram Mathematica таковы, что матрицы собственных векторов BE,n нужно транспонировать, прежде чем присвоить им имена WE,n.
Матрицы KEI и KEII, подобно Kx2, являются диагональными матрицами элементов, соответственно,
kz
n( -1 nI Sin
in (e)-i ij
ko
nil - I „i Sin
in (e)-i ЛГ
Нумерация индексов / в элементах К х , КЕд и КЕП проходит сверху вниз от /тш до ¿тах. Для четных 5: ¿тп = 1-5 / 2, /тах = 5-5 / 2; для нечетных 5: /тш = 1-(5+1) / 2, ¿тах = 5-(5+1) / 2.
При реализации ЦПМП, как показывалось выше, для каждого слоя (1 < п < N и области II (п = N+1) вводятся столбцы-константы Тп, через которые определяются столбцы-константы больших номеров или амплитуды прошедших дифракционных порядков, когда п = N. В матричном уравнении граничных условий на столбцы этих констант умножается вспомогательная матрица. Для п = N+1 она определяется из граничных условий, для остальных п - через еще одну вспомогательную матрицу, которая обеспечивает переход от слоя с большим номером к меньшему на 1. Таким образом через задание кусочной функции можно определить первую
I
jK E,II
We,„ (I + ХИЬИа-1ХИ)
VE,n ( — ХИМ-1Xn )
n = N +1
n Ф N +1
(17)
и вторую
an "We,„ Wen ' -1 fn+1
bn _ _ Ve n -Ve n _ _g n+1
(18)
вспомогательные матрицы. Поскольку (17) выражается через верхнюю an и нижнюю bn половины (18) отдельно, необходимо использовать функцию изъятия части списка. При этом обращение матрицы an (аич) задается отдельной функцией.
Далее следует цикл для n от N до n-1, в котором выполняется расчет вспомогательных матриц, их частей и инверсий an, поскольку матрицы для меньших номеров определяются через матрицы для больших номеров. После завершения цикла отдельно рассчитывается матрица (17) для n = 1, поскольку для первого слоя матрица (18) не требуется.
После определения f1 и g1 с помощью внутренних средств Wolfram Mathematica решается система линейных уравнений (10). Из матрицы-решения выводится нижняя половина-столбец (что соответствует T1), необходимая для расчета амплитуд прошедших дифракционных порядков по формуле (11), через операцию цикла. Для этого предварительно вводится функция
Ln
I
-'n+1Лn ^n
n = N +1 n Ф N +1
в цикле для п от N до 1, дающая множитель, состоящий из скалярных произведений перед Т в (11). Таким образом определяется матрица амплитуд прошедших дифракционных порядков Т.
Алгоритм завершается расчетом ДЭ необходимого порядка I по формуле:
DEt = TT* Re
V
nii - I nI sin (e)- i
nIcos0
(19)
v у
При алгоритмизации ГС-подхода используются граничные условия, содержащие диагональные матрицы K'I и К'п с элементами соответственно
. 2п kzIi . 2п 2 ( . . X j IT-f =j Tini —I "Is,n <e)—' Л
j T if = j TjЧ nisin (e)—i £)•
Отдельно задается функция от n
k
о
k
кг,п = £0,п — (1п(0))2.
X '
Матрица В'Е,п (формула (15)) строится следующим образом. Во-первых, создаются функции
0 а = 0
Еп = 2/кг,п;
во-вторых, соединяются четыре матрицы, причем нижняя левая Ап представляет собой сумму диагональной матрицы с элементами bi ( от ¿т1п до ¿тах ) и матрицы элементов апг_г, где г - номер строки, t - номер столбца. Правая нижняя Оп - диагональная матрица коэффициентов Еп, левая верхняя - нулевая матрица размером 5x5, правая верхняя I - единичная матрица.
Так же, как и в предыдущем подходе, создаются списки всех матриц В'Еп, их транспонированных матриц собственных векторов W'E,n и диагональных матриц их собственных чисел Q'n. Из последних двух отбираются элементы, соответствующие №прт1 и рпт, где 1 <р < 5 определяется через 5:
5
i + — для четных 5 2
Р = 1 1
5 +1
1 +--для нечетных 5.
I 2
Для построения итоговой матрицы граничных условий (15) создаются двумерные списки о с элементами С1,т,п = (Рп,т -/кг,п), 02 с элементами с2,т,п = ехр(с!,тп к-п / Щ, о3 с элементами с3,т,п = о2,т,п-ехр(с!,т,п к-п/ Щ. Далее с помощью них строятся матрицы ф1п, ф2,п, ф3,п, ф4,п и Фп в соответствии с выражением (15). В последней нулевые матрицы представляют собой двумерные списки нужной размерности со всеми нулевыми элементами. Элементы ф1п:
(-1) 1,п №пртЪз,т,п, элементы ф2,п:
( —1)81,п +1
элементы фз,п: ^1,т,п, элементы ф4,п:
С1,т,п С2,т,п, а Фп образуется через операцию соединения списков.
Из-за операции численного деления экспонент произведение с1т1 с3т1 в матрице ф21 может не дать единицу. Таким образом, ф2,1 задается отдельно как
( —1)81п +1 р т (Рт,п — /кг,п ) е'"" ^") ^ .
Матрицы 4, Z и v из (16) легко строятся соединением списков. Затем, также через соединение списков, но уже второго уровня, строится Н матрица. Средствами Wolfram Mathematica напрямую решается матричное уравнение (15) относительно всех неизвестных констант (ПРСУ-подход), и для определения амплитуд прошедших дифракционных порядков выбирают s нижних переменных из столбца-решения. То есть, поскольку вся система имеет 2s (N+1) переменных, выбираются элементы списка решений с номерами от s (2N+1)+1 до 2s (N+1). ДЭ рассчитывается по формуле (19).
Количество коэффициентов в (15) можно уменьшить с помощью Гауссовых сокращений (ГС-подход) [5]. Рассмотрим матрицы Фь 4, Z и v. По правилам матричных преобразований можно превратить вертикальное соединение ф^ и ф21 в матрицу с единицами на главной диагонали (тем самым -I и -K' в 4, 5i0 и -K'Se в v изменятся на (-I)(1) и (-K')(1) в 4 и (5Ю)(1) и (-K'5l0)(1) в v, где «(1)» символизирует изменение после исключения констант первого слоя). Затем следует обнулить все коэффициенты в соединении Ф3,1 и Ф4,1 (тем самым нулевые части в 4 и v, принадлежащие тем же строкам, что и коэффициенты Ф3,1 и Ф4,1, заполнятся сверху вниз матрицами (-I)(1), (-K')(1), (5i0)(1) и (-K'5i0)(1)). Отметим, что превращение в единичные матрицы лучше осуществлять без встроенной функции приведения к ступенчатому виду в системе Wolfram Mathematica, поскольку при больших s обнуляются нижние строки получившихся матриц.
Таким образом, все константы для первого периодического слоя не примут участия в дальнейших расчетах. В получившемся измененном матричном уравнении (15) можно повторить эту операцию для остальных периодических слоев. В результате останется матричное уравнение, содержащее только амплитуды отраженных и прошедших дифракционных порядков, то есть
(-I)(N) -I ■ " R" ■ (-0 Г "
(-KI ) > KIi _ T (-KI 5,o ) >_
Отбор нижней половины столбца-решения даст амплитуды прошедших дифракционных порядков. Дифракционная эффективность рассчитывается по формуле (19).
5. Результаты и обсуждение Для анализа выбиралась структура с параметрами: п = 1; пп = 1,5; высота рельефа к = 1 мкм; параметр Р = Л / к = 5; число слоев N = 40; случай ТЕ-поляризации. Сравнивалась ДЭ в 1-м порядке дифракции при числе учитываемых дифракционных порядков 5 = 4, 8, 20, 40. В табл. 1 приведены результаты, полученные при использовании для расчета ДЭ трех подходов: ЦПМП-подхода, прямого решения системы уравнений граничных условий (ПРСУ) и ГС-подхода.
Табл. 1. Результаты расчета ДЭ, полученные с использованием ЦПМП-подхода, путем прямого решения системы уравнений граничных условий (ПРСУ) и ГС-подхода при различных значениях параметров: Г с нижними индексами - время, затраченное на расчет с использованием соответствующего подхода
« = 1; nII = 1,5; h = 1 мкм; P = Л / h = 5; N = 40; i = 1
s = 4, X = 300 нм s = 20, X = 300 нм
ö,° ЦПМП ПРСУ ГС ЦПМП ПРСУ ГС
ДЭ ДЭ ¿прсу/ ¿цпмп ДЭ ¿гс/ ¿цпмп ДЭ ДЭ ¿прсу/ ¿цпмп ДЭ ¿гс/ ¿цпмп
-30 0,386786 0,386786 3 0,386786 7,6 0,197272 0,197203 6,1 0,197554 21,17
-15 0,468206 0,468206 8,1 0,468206 9,1 0,232136 0,232136 6,89 0,232136 18,17
0 0,49717 0,49717 5,1 0,49717 5,8 0,196175 0,196175 6,85 0,194886 22,11
15 0,459018 0,459018 7,2 0,459018 11 0,101957 0,101957 7,29 0,101965 22,46
30 0,340033 0,340033 6,8 0,340033 7,8 0,011972 0,0117936 7,81 0,013911 20,27
s = 4, X = 400 нм s = 20, X = 400 нм
ö,° ЦПМП ПРСУ ГС ЦПМП ПРСУ ГС
ДЭ ДЭ ¿прсу/ ¿цпмп ДЭ ¿гс/ ¿цпмп ДЭ ДЭ ¿прсу/ ¿цпмп ДЭ ¿гс/ ¿цпмп
-30 0,599701 0,599701 6,1 0,599701 9,6 0,680996 0,685046 7,11 0,685727 20,68
-15 0,676283 0,676283 5,7 0,676283 8 0,760327 0,760327 6,77 0,760213 20,34
0 0,722416 0,722416 5,9 0,722416 9,3 0,753267 0,753267 6,44 0,753045 21,33
15 0,737819 0,737819 7,2 0,737819 7.1 0,648659 0,648632 6,94 0,648683 20,63
30 0,690355 0,690355 6,1 0,690355 5,7 0,425972 0,434048 6,34 0,393695 18,5
s = 4, X = 500 нм s = 20, X = 500 нм
ö,° ЦПМП ПРСУ ГС ЦПМП ПРСУ ГС
ДЭ ДЭ ¿прсу/ ¿цпмп ДЭ ¿гс/ ¿цпмп ДЭ ДЭ ¿прсу/ ¿цпмп ДЭ ¿гс/ ¿цпмп
-30 0,613634 0,613634 5,9 0,613634 9,1 0,791659 0,843981 6,55 0,851353 24,81
-15 0,656357 0,656357 5,7 0,656357 8,4 0,861832 0,863335 7,19 0,863333 20,9
0 0,696179 0,696179 5,5 0,696179 8,25 0,880109 0,880109 6,87 0,880132 18,85
15 0,740601 0,740601 7,6 0,740601 9,6 0,853823 0,853607 5,82 0,851754 22,1
30 0,770812 0,770812 5,4 0,770812 6,4 0,739961 0,74614 7,17 1750,18 21,67
s = 8, X = 300 нм s = 40, X = 300 нм
ö,° ЦПМП ПРСУ ГС ЦПМП ПРСУ ГС
ДЭ ДЭ ¿прсу/ ¿цпмп ДЭ ¿гс/ ¿цпмп ДЭ ДЭ ¿прсу/ ¿цпмп ДЭ ¿гс/ ¿цпмп
-30 0,230078 0,230078 5,4 0,230078 14,25 0,1898 0,17463 8,66 4,3-10" 36,88
-15 0,280596 0,280596 5,6 0,280596 14,75 0,22422 0,221593 8,97 1294,71 36,76
0 0,258874 0,258874 5,2 0,258874 11,2 0,190739 0,1898 8,49 0,16973 35,97
15 0,15461 0,15461 5,8 0,15461 14,75 0,097263 0,0976295 9,14 14,2653 39,32
30 0,0277394 0,02773 4,8 6,19236 14,25 0,010753 0,0102328 8,59 2,81015 38,02
s = 8, X = 400 нм s = 40, X = 400 нм
ö,° ЦПМП ПРСУ ГС ЦПМП ПРСУ ГС
ДЭ ДЭ ¿прсу/ ¿цпмп ДЭ ¿гс/ ¿цпмп ДЭ ДЭ ¿прсу/ ¿цпмп ДЭ ¿гс/ ¿цпмп
-30 0,679211 0,679211 4,33 0,679211 11,2 0,68121 0,718119 8,92 2,31012 37,66
-15 0,762063 0,762063 4,5 0,762063 11,6 0,755157 0,779014 9,13 1,5-109 37,73
0 0,778282 0,778282 5 0,850481 10,8 0,748268 0,75885 8,64 265,822 38,52
15 0,706478 0,706478 5,6 0,706478 14,75 0,644913 0,651165 8,59 2-10° 38,25
30 0,490053 0,490053 5,2 0,490053 13,75 0,420677 0,433868 8,93 4,51020 37,74
s = 8, X = 500 нм s = 40, X = 500 нм
ö,° ЦПМП ПРСУ ГС ЦПМП ПРСУ ГС
ДЭ ДЭ ¿прсу/ ¿цпмп ДЭ ¿гс/ ¿цпмп ДЭ ДЭ ¿прсу/ ¿цпмп ДЭ ¿гс/ ¿цпмп
-30 0,753371 0,753371 4,5 0,753371 11,6 0,796477 0,891731 8,28 1,3-10" 37,02
-15 0,816468 0,816468 4,7 0,816468 12 0,858472 0,895091 9,17 41014 38,92
0 0,846104 0,846104 5,2 0,846104 14,25 0,876235 0,890746 8,37 109 37,82
15 0,832215 0,832215 5,8 0,832215 10,17 0,853523 0,853013 9,12 1013 37,16
30 0,742428 0,742428 11,6 0,742424 11,4 0,741707 0,746174 7,94 1,5-10" 37,3
Анализ совпадения или несовпадения значений ДЭ в рамках разных подходов, но при одинаковых параметрах X и 0 производился для ряда различных значений длин волн (X = 300 нм, 400 нм, 500 нм) и углов падения излучения на микроструктуру (0 = -30°, -15°, 0°, 15°, 30°). Для анализа производительности RCWA-метода, базирующегося на том или ином подходе, в табл. 1 приведены средние затраты времени на расчет ДЭ, нормированные на затраты времени, требовавшиеся при расчете, базировавшемся на ЦПМП-подходе.
Как видно из табл. 1, значения ДЭ, полученные с помощью алгоритмов, основанных на ЦПМП- и ПРСУ-подходах, идентичны для s = 4 и s = 8. Различия появляются при s = 20 и s = 40, что объясняется плохой обусловленностью (почти вырожденностью) больших матриц коэффициентов (см. выражение (15)) при неизвестных в ПРСУ. Отметим, что среда Wolfram Mathematica зафиксировала эту плохую обусловленность. Чем хуже обусловлена матрица коэффициентов при решении системы линейных алгебраиче-
ских уравнений, тем выше расхождения найденных и истинных значений искомых величин.
Такая же проблема имеет место и при использовании ГС-подхода. Действительно, после увеличения параметра 5 до определенного значения или при определенных параметрах X и 0 и небольшом 5 алгоритм становится неустойчивым, о чем говорят значения ДЭ, превышающие единицу, что нефизично (см. табл. 2). Значения 5, приводящие к неустойчивости, зависят от параметров ДОЭ. Как показывает табл. 2, для 0 = 0° и X = 400 нм неустойчивость резко возникает при 5 = 33, а для 0 = 30° и X = 300 нм она локально появляется при 5 = 8, в то время как при других близких значениях 5 ДЭ сравнима с полученной с использованием подходов ЦПМП и ПРСУ (см. табл. 1). Табл. 2. Достижение неустойчивости алгоритмом ГС при возрастании числа порядков дифракции 5 в анализе
nI = 1; nII = 1,5; h = 1 мкм; P = Л / h = 5; N = 40; i = 1; 0 = 0°; X = 400 нм
s 25 26 27 28
ДЭ 0,7507 0,7506 0,749 0,7487
s 29 30 31 32
ДЭ 0,759 0,745 0,75 0,747
s 33 34 35 36
ДЭ 4,325 0,4856 5,5723 33,29
s 37 38 39 40
ДЭ 147,7 249,3 4260,13 265,822
nI = 1; nII = 1,5; h = 1 мкм; P = Л / h = 5; N = 40; i = 1; 0 = 30°; X = 300 нм
s 3 4 5 6
ДЭ 0,4405 0,34 0,0722 0,0585
s 7 8 9 10
ДЭ 0,0221 6,1924 0,0211 0,0205
s 11 12 13 14
ДЭ 0,0166 0,0197 0,0117 0,0136
s 15 16 17 18
ДЭ 0,0064 0,0108 0,0123 0,0123
Неустойчивость в таких случаях обусловлена неудавшимся приведением соединений матриц Ф1,1 и Ф2,1 к матрицам с единицами на главной диагонали. Это происходит из-за недостаточной точности вычислений (большой величины машинного нуля) в системе Wolfram Mathematica и в большинстве других сред - 10-16. Поскольку матрицы Фя тем больше по размеру, чем выше s, то после операций сокращения последовательно для каждого слоя обратные элементы в измененных 4 и v все ближе к машинному нулю. Соответственно, точно выполнить необходимые вычитания для приведения матриц Фя к нужному для исключения неизвестных констант виду не представляется возможным численно.
Возможен вариант, когда используется округление коэффициентов в Фя, 4 и v с заданной точностью (Wolfram Mathematica способна округлять и приводить числа к обыкновенной дроби), что позволяет обойти ограничение максимальной точности чисел с плавающей запятой. Однако при этом значительно увеличивается время расчета, а точность вычислений и устойчивость могут снизиться из-за численных усечений, обусловленных ошибками округления.
Сравнение временных интервалов, затрачиваемых на расчет ДЭ, показывает, что наиболее эффективен в этом отношении алгоритм, базирующийся на ЦПМП-подходе, благодаря использованию оптимизированных внутренних алгоритмов обращения матриц. Время, затрачиваемое на расчет ДЭ, при использовании ПРСУ-подхода больше из-за размера системы уравнений и ее непосредственного решения. При реализации ГС-подхода время на расчет по сравнению с алгоритмом, базирующимся на ПРСУ-подходе, увеличивается с ростом числа учитываемых дифракционных порядков по причине возрастания размера системы уравнений, подвергающейся многочисленным матричным преобразованиям.
Результаты, представленные в табл. 1, показывают, что при увеличении числа учитываемых дифракционных порядков ДЭ асимптотически приближается к определенному значению, которое можно назвать «истинным» в рамках данного метода. Например, алгоритмы, базирующиеся как на ЦПМП-, так и на ПРСУ-подходе, дают при 5 = 4 и 5 = 8 существенно различающиеся значения ДЭ, в то время как при 5 = 20 и 5 = 40 разница в значениях ДЭ не превышает 0,07.
Приведенный в данном параграфе анализ показывает, что ГС-подход является наименее устойчивым и не имеющим значимых преимуществ по сравнению с другими подходами. Поэтому его нецелесообразно использовать для вычисления ДЭ при большом количестве учитываемых дифракционных порядков и значительном числе слоев в микроструктуре. ПРСУ-подход более стабилен, чем ГС-подход, но менее стабилен и эффективен по сравнению с ЦПМП-подходом.
В результате наиболее подходящим алгоритмом для решения задач, предполагающих расчет ДЭ дифракционной микроструктуры, глубина рельефа которой одного порядка с пространственным периодом, а число слоев в рельефе относительно не велико, но оптимально, является алгоритм, основанный на ЦПМП-подходе.
Литература
1. Грейсух, Г.И. Минимизация суммарной глубины внутренних пилообразных рельефов двухслойной рельефно-фазовой дифракционной микроструктуры / Г.И. Грейсух, В.А. Данилов, С.А. Степанов,
A.И. Антонов, Б.А. Усиевич // Оптика и спектроскопия. - 2018. - Т. 124, Вып. 1. - С. 100-104.
2. Грейсух, Г.И. Спектральная и угловая зависимости эффективности трехслойных рельефно-фазовых дифракционных элементов ИК диапазона / Г.И. Грейсух, В.А. Данилов, С.А. Степанов, А.И. Антонов, Б.А. Усиевич // Оптика и спектроскопия. - 2018. - Т. 125, Вып. 1. - С. 57-61.
3. Грейсух, Г.И. Гармоническая киноформная линза: дифракционная эффективность и хроматизм / Г.И. Грейсух,
B.А. Данилов, С.А. Степанов, А.И. Антонов, Б.А. Усиевич // Оптика и спектроскопия. - 2018. - Т. 125, Вып. 2. -
C. 223-228.
4. Грейсух, Г.И. Спектральная и угловая зависимость эффективности двухслойной однорельефной пилообраз-
ной микроструктуры / Г.И. Грейсух, В.А. Данилов, А.И. Антонов, С.А. Степанов, Б.А. Усиевич // Компьютерная оптика. - 2018. - Т. 42, № 1. - С. 38-43. - DOI: 10.18287/2412-6179-2018-42-1-38-43.
5. Moharam, M.G. Diffraction analysis of dielectric surface-relief gratings / M.G. Moharam, T.K. Gaylord // Journal of the Optical Society of America. - 1982. -Vol. 72, Issue 10. - P. 1385-1392.
6. Moharam, M.G. Formulation for stable and efficient implementation ofthe rigorous coupled-wave analysis of binary gratings / M.G. Moharam, E.B. Grann, D.A. Pommet, T.K. Gaylord // Journal of the Optical Society of America. -1995. - Vol. 12, Issue 5. - P. 1068-1076.
7. Moharam, M.G. Stable implementation of the rigorous coupled-wave analysis for surface-relief gratings: Enhanced transmittance matrix approach / M.G. Moharam, E.B. Grann, D.A. Pommet, T.K. Gaylord // Journal of the Optical Society of America. - 1995. - Vol. 12, Issue 5. - P. 1077-1086.
8. Li, L. Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings / L. Li // Journal of the Optical Society of America. - 1996. -Vol. 13, Issue 6. - P. 1024-1035.
9. Redheffer, R. Difference equations and functional equations in transmission-line theory // R. Redheffer / Modern Mathematics for the Engineer. - 1961. - Vol. 12. - P. 282-337.
10. CEM Lectures [Электронный ресурс]. - 2018. - URL: http://emlab.utep.edu/academics.htm (дата обращения 1.09.2018).
11. Солимено, С. Дифракция и волноводное распространение оптического излучения : пер. с англ. / С. Солимено, Б. Крозиньяни, П. Ди Порто. - М.: Мир, 1989. - 664 с.
12. Тихонов, А.Н. Дифференциальные уравнения /
A.Н. Тихонов, А.Б. Васильева, А.Г. Свешников. - М.: Физматлит, 2005. - 256 с.
13. Кирьянов, Д.В. Вычислительная физика / Д.В. Кирьянов, Е.Н. Кирьянова. - М.: Полибук Мультимедиа, 2006. - 352 с.
14. Ильин, В.А. Основы математического анализа: в 2-х ч. /
B.А. Ильин. - Часть II. - М.: Физматлит, 2002. - 464 с.
15. Федеральное государственное бюджетное учреждение науки. Институт общей физики им. А.М. Прохорова Российской академии наук [Электронный ресурс]. -2018. - URL: www.gpi.ru (дата обращения 1.09.2018).
16. Grating Diffraction Calculator [Электронный ресурс]. -2018. - URL: kjinnovation.com/GD-Calc.html (дата обращения 1.09.2018).
17. MATLAB [Электронный ресурс]. - 2018. - URL: uk.mathworks.com/products/matlab.html (дата обращения 1.09.2018).
18. Wolfram Mathematica: Современные технические вычисления [Электронный ресурс]. - 2018. - URL: www.wolfram.com/mathematica.html (дата обращения 1.09.2018).
Сведения об авторах
Антонов Артем Иванович, 1993 года рождения. В 2014 окончил бакалавриат Пензенского государственного университета по направлению «Физика». В 2016 окончил магистратуру Пензенского государственного университета по направлению «Физика». Аспирант кафедры физики и химии Пензенского государственного университета архитектуры и строительства по направлению «Приборы и методы экспериментальной физики». В списке научных работ А.И. Антонова 4 статьи и 1 авторское свидетельство. E-mail: [email protected] .
Васин Леонид Анатольевич. 1971 года рождения. В 1996 году окончил Пензенский технологический институт (завод-ВТУЗ) по специальности «Вычислительные машины, системы и сети». Кандидат технических наук (2002 год). Работает заведующим кафедрой информационно-вычислительных систем. В списке научных трудов более 40 статей, один патент. E-mail: [email protected] .
Грейсух Григорий Исаевич, 1943 года рождения. В 1965 году окончил Пензенский политехнический институт по специальности «Радиотехника». Заслуженный работник высшей школы РФ, доктор технических наук (1990 год), профессор. Работает заведующим кафедрой физики и химии Пензенского государственного университета архитектуры и строительства. Является членом Европейского оптического общества (EOS) и Российского оптического общества им. Д.С. Рождественского. Г.И. Грейсух - специалист в области расчета оптических систем, дифракционной и градиентной оптики. В списке научных работ Г.И. Грейсуха более 200 статей, 3 монографии, 9 авторских свидетельств. E-mail: [email protected] .
ГРНТИ: 29.31.15
Поступила в редакцию 25 января 2019 г. Окончательный вариант - 16 марта 2019 г.
Approaches to the algorithmization of the rigorous coupled-wave analysis
A.I. Antonov1, L.A. Vasin1, G.I. Greisukh1 1Penza State University of Architecture and Construction, Penza, Russia
Abstract
Coefficients required for the algorithmization and implementation of the rigorous coupled-waves analysis are derived for a Fourier series expansion of the permittivity of a linear-sawtooth relief with positive tangent of the working surface angle and vertical backward slope. Approaches to the implementation of the rigorous coupled-waves analysis are described and compared in terms
of stability and efficiency, namely, the approach of a purposeful transformation of the transmission matrix and the approach of Gaussian illuminations. It was concluded that the most appropriate approach for determining the diffraction efficiency for cases of deep diffraction structures is the approach of purposeful transformation of the transmission matrix.
Keywords: diffraction efficiency, permittivity, Maxwell's equations, rigorous coupled-wave analysis.
Citation: Antonov AI, Vasin LA, Greisukh GI. Approaches to the algorithmization of the rigorous coupled-wave analysis. Computer Optics 2019; 43(2): 209-219. DOI: 10.18287/2412-61792019-43-2-209-219.
References
[1] Greisukh GI, Danilov VA, Stepanov SA, Antonov AI, Usievich BA. Minimization of the total depth of internal sawtooth reliefs of a two-layer relief-phase diffraction microstructure. Optics and Spectroscopy 2018; 124(1): 98-102.
[2] Greisukh GI, Danilov VA, Stepanov SA, Antonov AI, Usievich BA. Spectral and angular dependences of efficiency of three-layer relief-phase diffraction elements of the IR range. Optics and spectroscopy 2018; 125(1): 60-64.
[3] Greisukh GI, Danilov VA, Stepanov SA, Antonov AI, Usievich BA. Harmonic kinoform lens: Diffraction efficiency and chromatism. Optics and spectroscopy 2018; 125(2): 232-237.
[4] Greisukh GI, Danilov VA, Stepanov SA, Antonov AI, Usievich BA. Spectral and angular dependence of the efficiency of a two-layer and single-relief sawtooth microstructure. Computer Optics 2018; 42(1): 38-43. DOI: 10.18287/2412-6179-2018-42-1-38-43.
[5] Moharam MG, Gaylord TK. Diffraction analysis of dielectric surface-relief gratings. J Opt Soc Am 1982; 72(10): 1385-1392.
[6] Moharam MG, Grann EB, Pommet DA, Gaylord TK. Formulation for stable and efficient implementation ofthe rigorous coupled-wave analysis of binary gratings. J Opt Soc Am 1995; 12(5): 1068-1076.
[7] Moharam MG, Grann EB, Pommet DA, Gaylord TK. Stable implementation of the rigorous coupled-wave analysis
for surface-relief gratings: enhanced transmittance matrix approach. J Opt Soc Am 1995; 12(5): 1077-1086.
[8] Li L. Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings. J Opt Soc Am 1996; 13(6): 1024-1035.
[9] Redheffer R. Difference equations and functional equations in transmission-line theory. Modern Mathematics for the Engineer 1961; 12: 282-337.
[10] CEM Lectures. Source: (http://emlab.utep.edu/academics.htm).
[11] Solimeno S, Crosignani B, DiPorto P. Guiding, diffraction, and confinement of optical radianion. Orlando, FL: Academic Press Inc; 1989.
[12] Tikhonov AN, Vasil'eva AB, Sveshnikov AG. Differential equations. Berlin, Heidelberg: Springer-Verlag; 1985.
[13] Kiryanov D, Kiryanova E. Computational science. Hung-ham, MA: Infinity Science Press LLC; 2007.
[14] Ilyin VA, Poznyak EG. Fundamentals of mathematical analysis: Volumes 1 and 2. Part II. Moscow: "Fizmatlit" Publisher; 1982.
[15] Prokhorov General Physics Institute. Source: (http://www.gpi.ru).
[16] Grating Diffraction Calculator. Source: (http://kjinnovation.com/GD-Calc.html).
[17] MATLAB. Source: (https://uk.mathworks.com/products /matlab.html).
[18] Wolfram Mathematica. Source: (http://www.wolfram.com /mathematica.html).
Author's information
Artyom Ivanovich Antonov (b. 1993). In 2014 he graduated from Penza State University, direction of "Physics", as bachelor. In 2016 he acquired the master degree of direction of "Physics". Now he is a post-graduate student of Physics and Chemistry department of the Penza State University of Architecture and Construction in the direction "Devices and Methods of Experimental Physics". A.I. Antonov is co-author of 4 articles and 1 invention. E-mail: [email protected] .
Leonid Anatolyevich Vasin (b. 1971). In 1996, he graduated from Penza Technological Institute (plant-VTUZ) with a degree in Computing Machines, Systems and Networks. Candidate of Technical Sciences (2002). Works as head of Information and Computing Systems department. In the list of scientific works of more than 40 articles, one patent. E-mail: [email protected] .
Grigoriy Isaevitch Greisukh (b. 1943) graduated (1965) from Penza Politechnical Institute, majoring in Radio Engineering. He is the deserved worker of the Russian Higher School. He received his Doctor in Technical (1990) degrees from the Leningrad Institute of Precision Mechanics and Optics. He is chief of the Physics and Chemistry department of the Penza State University of Architecture and Construction. G.I. Greisukh is EOS and D. S. Rozhdestvensky Optical Society member. His current research interests include design of optical system, diffractive and gradient-index optics. He is co-author of 200 scientific papers, 3 monographs, and 9 inventions. E-mail: [email protected] .
Received January 25, 2019. The final version - March 16, 2019.
