Компьютерный расчет рассеяния электромагнитных волн на одномерных дифракционных решетках методом матричного уравнения Риккати

«Московский физико-технический институт (национальный исследовательский университет)», 2019

связь, оптические вычисления, биосенсорика и т.д., в том числе при создании высокопроизводительных фотонных аналого-цифровых преобразователей [2]. Базовыми компонентами любого радиофотонного устройства являются волноводы, служащие как для передачи света, так и для его преобразования и обработки. Соответственно, возникает задача ввода электромагнитного (ЭМ) излучения в волновод или вывода из него. Среди методов решения этой задачи весьма перспективным является дифракционный метод, предполагающий использование различных периодических структур в составе самого волновода. Вследствие малых и сравнимых с длиной волны вводимого/выводимого излучения размеров рассматриваемых устройств необходим учет различных явлений взаимного преобразования однородных (распространяющихся) и неоднородных (затухающих) волн, возникающих при их многократном рассеянии на неоднородностях границ раздела и объема сред. Таким образом, становится актуальной задача рассеяния света на таких неоднородных структурах.

Исторически сложилось, что в теории многократного рассеяния волновых полей независимо развивались методы описания рассеяния на неоднородных, в том числе периодических, поверхностях и в неоднородных объемных диэлектрических средах. Рассеяние на поверхности обычно описывается в терминах граничной задачи связанной в общем случае с векторным волновым уравнением, но в некоторых случаях сводящейся к скалярному уравнению Гельмгольца с соответствующими граничными условиями. Методы решения этой задачи основываются по большей части на интегральных уравнениях, записываемых для граничного поля, с последующим их численным решением [3]. Также применяется техника матрицы рассеяния (S-матрицы). В теории многократного рассеяния в объемных средах известны три подхода: метод композиции операторов рассеяния Ватсона (метод Т-мат-рицы) [3], метод инвариантного погружения [4] и метод трансфер-матриц (матриц переноса) [5].

Исследования фотонных кристаллов привели к возникновению новых способов расчета коэффициентов отражения и прохождения волн для периодических объемных структур, вычисления локальных полей в них и расчета зонного характера их частотного спектра, появление которого обусловлено геометрическим резонансом при брэгговском рассеянии волн. Таким методом является метод плоских волн (см., например, [6]), который сводится к матричному интегральному уравнению Фредгольма, что характерно для задач поверхностного рассеяния в их традиционной формулировке, а также метод конечных разностей в различных вариациях [3, 7]. Последний вместе с техникой S-матрицы близок к методу инвариантного погружения, и потому эти подходы могут быть рассмотрены как шаг к единообразному описанию поверхностного и объемного рассеяния ЭМ волн. Следующим решающим шагом явились так называемые «соотношения переноса» [8] в теории многократного рассеяния волн, которые связывают между собой локальные поля внутри слоя рассеивающей среды и коэффициенты отражения и прохождения этого слоя. На основе соотношений переноса развит единый точный подход [9] к теоретическому описанию многократного резонансного рассеяния волн на одномерных периодических поверхностях раздела двух диэлектрических сред (оптических решетках) и на двумерных периодических структурах. Основным инструментом упомянутого точного подхода, дающего возможность численного описания свойств резонансных решеток и фотонных структур, является матричное уравнение Риккати для коэффициента отражения и ассоциированное с ним уравнение для коэффициента прохождения периодической структуры. Для краткости будем называть этот подход методом уравнения Риккати. Важной особенностью является то, что эти уравнения выведены с учетом взаимного преобразования однородных и неоднородных волн, и, таким образом, этот метод позволяет преодолевать трудности, связанные с геометрическим расположением рабочей области микрофотонных устройств в ближней зоне дифракционной решетки.

• 2.    Метод матричного уравнения Риккати

Рассмотрим рассеяние плоской монохроматической волны на абстрактной неоднородной периодической структуре (рис. 1), локализованной в слое толщиной Һ в декартовой системе координат и характеризуемой также периодом Л и диэлектрической проницаемостью г(х,г,ш ), г де ш — частота ЭМ излучения. Плоская линейно поляризованная волна с длиной волны в вакууме А падает из верхнего полупространства под углом 9 к оси г. Рассматривается случай поперечно-электрической (ТЕ) поляризации, что позволяет свести векторное волновое уравнение к скалярному уравнению Гельмгольца для комплексных амплитуд.

Рис. 1. Иллюстрация к описанию метода уравнения Риккати

Метод позволяет представить поле отраженного и прошедшего излучений в виде совокупности волн с определенными волновыми векторами k± = (^_м, 0, ±а ^ )^т ^ компоненты которых определяются параметрами исходного излучения и структуры:

2^^

Р ц = ко sin 9 + —,

д = 0, ±1, ±2,...,

а? = У ^{к 2} - $    Д = 0, ±1, ±2,...,                      (2)

где к = 2^л/Ё/А — волновое число в среде с диэлектрической проницаемостью s, в частности к₀ — волновое число в св ободном пространстве (s = Еь_д = 1), обозначение а_ц без верхнего индекса соответствует ко. Волну, имеющую волновой вектор k±, будем называть ц-й модой решетки в спектрах отражения и прохождения соответственно, при этом нулевая мода в спектре отражения определяет зеркальное направление. Волновые векторы в общем случае комплексные, так как комплексной может быть как диэлектрическая проницаемость среды, так и одна из компонент вектора, а именно г-компонента. Причем, если выражение под корнем станет отрицательным, величина а окажется чисто мнимой, а это значит, что соответствующая волна перейдет из распространяющейся (однородной) в затухающую (эванесцентную, неоднородную). В таком случае модуль а называют коэффициентом затухания.

Уравнение Риккати, записанное для матричного коэффициента отражения Д, вместе с граничным (или «начальным») условием имеет следующий вид:

— = ЯА(г)Д + (Л(г) + С )Д + Я(А(г) + С) + А(г),    0 < г 6 Һ,

R|,^o = Ro, где R, Ro, А, С — бесконечномерные квадратные матрицы, компоненты которых нумеруются индексами ц, v = 0, ±1, ±2,...

Граничным (начальным) условием является матричный коэффициент отражения от однородного слоя подложки, определяющийся выражением

(W = (Rn/)ц

1 - exp (2г<71ц Lo)    ^

1 - (Rn/)Ц exp (2iai_ll Lo)

^v = 0, ±1, ±2,...,

где Фр. ~ 2-компонента ц-го волнового вектора в подложке, Lo — толщина подложки, \ 0, ц = v o^v = <          — символ Кронекера, a (Rjn/)ц представляет собой коэффициент отра-

• 1,    ц = ^v

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

£ц_-_£1ц

^(R inf )Ц          .       .

СТ р + ^щ

Входящая в уравнение (3) диагональная матрица С образована 2-компонентами волновых векторов k⁺ излучения, рассеянного в v-й порядок, домноженными на мнимую единицу:

^С рр — ⁱ^ p d ₉ , ,     ^ц, ^v — ⁰, ±1, ±^2,...

Матрица А(2) включает в себя потенциал рассеяния переходной области, а ее компоненты определяются выражениями кп

^А рр ⁽²) = ⁱ д-^р - р (²),      ц^,v = 0, ±1, ±^2,...,

2^^ц где функция /^(2) носит название функции трансформации и описывает взаимное влияние однородных и эванесцентных волн при их многократном рассеянии в неоднородной области с периодом Л, и задается следующим образом:

Л

УД2) = Л / е(ж,2)е-2^ 0

dx.

Здесь интеграл берется внутри одного периода, зависимость диэлектрической проницаемости е от частоты ш ЭМ излучения опущена, так как последняя полагается фиксированной. В большинстве случаев простых симметричных структур интеграл удается легко вычислить аналитически, когда диэлектрическая проницаемость представляется кусочнооднородной функцией координат. Так, например, для дифракционной решетки с профилем сечения в виде равностороннего треугольника с высотой һ и основанием Л, с диэлектрической проницаемостью е в среде с е^д (рис. 2) функция трансформации принимает вид е — еЬд               2

УД2) =

— sin ^ц (1 — һ)).

Электрическое поле в полупространстве 2 > һ записывается в виде суммы полей падающей и отраженных волн:

∞

S(r) = exp (iko(r — hnz)) + ^ R^o(h)exp (ik+(r — hnz)) ,             (5)

ц=-^

где E(r) = (0,S(r), 0)T для ТЕ по.тяризащш. nz = (0, 0,1)т. r = (x,y,2)T. Таким образом, видим, что компоненты R^o(h) имеют смысл парциальных коэффициентов отражения падающей волны от периодической неоднородной структуры в ц-ю моду, образующуюся в неоднородной области при рассеянии на ней падающего излучения.

Аналогичным образом записывается ассоциированное с уравнением Риккати дифференциальное уравнение для матричного волнового коэффициента прохождения Т через неоднородную периодическую структуру и граничное (начальное) условие:

— = Т ⁽A(z) + С + A(z)R(z )⁾,    0 6 h,

Т |г ,: = То, где То представляет собой матричный коэффициент прохождения однородного слоя (подложки) толщиной L0, соответствующий матричному коэффициенту отражения R0 (4) от него, с компонентами

^(Т0 ) р.и ^—

_____________________4<Тл£1^_

(op + ощ)² exp (-гстщ L_o ) - (o ^

- аід)² exp (гоі л Lo )

S ^v

^,v — 0, ±1, ±2,...

Электрическое поле в полупространстве z 6 -Lo определяется выражением

∞

Е ( r ) — У До(h)exp(i k -( r + Lo n _z)).

л=—^

• 3.    Результаты расчета

• 3.1.    Серебряная решетка с треугольным профилем

  

  Рис. 2. Серебряная решетка с треугольным профилем

  
  

Для решения уравнений (3), (6) с соответствующими граничными условиями ограничим размерность входящих в них матриц некоторым числом М — 2п + 1, которое также будет иметь смысл количества мод. Это представляется возможным вследствие того, что моды высоких порядков имеют малые амплитуды и быстро затухают вдоль оси z, и, следовательно, дают малый вклад в суммы (5), (7). В настоящей работе производится учет 21 моды. Решение матричного уравнения Риккати для коэффициента отражения R осуществляется с помощью численных методов класса Рунге—Кутты, а именно вложенных методов 5-го порядка с оценкой локальной погрешности с помощью выражений 4-го порядка (Dormand— Prince, Cash—Karp). Для решения уравнения относительно коэффициента прохождения Т используется модифицированный метод Эйлера типа предиктор-корректор 2-го порядка. Рассматривается серебряная решетка с треугольным профилем сечения (рис. 2) и прямоугольная решетка в структуре кремний-на-изоляторе (рис. 5), расположенные на подложках (бесконечно) большой толщины.

Параметры излучения: длина волны А — 632.8 нм, угол падения д.

Параметры структуры: высота решетки h, пер иод Л, диэлектрическая проницаемость s( А) — -17.5 + 0.7г.

Рис. 3. Зависимости модулей коэффициентов отражения 1Н_ц0 | от отношения Х/Л. Вертикальными штриховыми линиями показаны моменты «открытия» мод — слева направо на левом рисунке ц = ±1, ±2, ±3, ±4. на правом рисунке ц = 1, 2, 3,4 — линии с коротким штрихом, ц = — 1, -2, —3 — с длинным. Стрелками указано, какая именно мода открывается. Сплошной и полой стрелками на правом рисунке обозначено смещение моментов открытия положительных и отрицательных порядков при изменении угла падения. Отрицательные порядки не представлены сами, однако оказывают воздействие на другие моды

При фиксированных длине волны излучения и высоте решетки ( Һ = 696.1 нм) получены зависимости коэффициентов отражения от отношения Х/Л при двух значениях угла падения р = 0, 0.4 рад (рис. 3).

. Х sinp +лц = 1,

ц = ц/іж = ±1, ±2,...

Резонансы, отображенные на рис. 3 и связанные с открытием мод, носят название аномалий Вуда (Рэлея—Вуда) [10].

Существует также резонансное поведение в зависимости от отношения высоты решетки к длине волны падающего излучения. Такие аномалии получили название параллельных аномалий Вуда или аномалий Вуда—Палмера [11]. Они продемонстрированы на рис. 4. Зависимости коэффициентов отражения получены при фиксированных длине волны излучения, угле падения ( р = 0) и периоде решетки (Л = 1 мкм).

Рис. 4. Зависимости модулей коэффициентов отражения \Нр,₀ | от высотьi решетки К. Синими прямоугольниками обозначена оценка Хесселя—Олинера для минимумов по отражению

• 3.2.    Кремниевая прямоугольная решетка на оксидном слое

Рис. 5. Кремниевая прямоугольная решетка с равномерным заполнением на оксидном слое

Параметры излучения: длина волны А = 1.55 мкм, угол падения ф = 0.

Параметры структуры: высота решетки һ, пер иод Л, фактор заполнения — отношение ширины «зуба» к периоду — Ғ = 0.5, толщина оксидного слоя L_OXl диэлектрические проницаемости кремния ^(А) = 12.041 и ок сида е _ох (А) = 2.085.

Параметры с указанными значениями считаются фиксированными.

При заданной высоте ( һ = 50 нм) построены зависимости коэффициентов отражения от отношения А/Л при двух значениях толщины оксидного слоя L_OX = 0, 0.2 мкм (рис. 6).

Рис. 6. Зависимости модулей коэффициентов отражения |Дмо | от отношения А/Л. Штриховыми линиями показаны моменты открытия мод — слева направо д = ±1, ±2, ±3, ±4, стрелками обозначены непосредственно открывающиеся моды (собственные резонансы)

На графиках также проявляются резонансы, связанные с открытием мод, но, кроме того, характер зависимостей значительно изменяется при добавлении слоя оксида кремния. Собственные резонансы становятся резче, как и их влияние на моды других порядков (см., например, открытие ±1 мод). Приведем далее распределения отраженного поля вдали от резонанса, когда однородной является только нулевая мода, (А/Л to 1.25) и непосредственно вблизи открытия ±1 мод (А/Л to 1.0) в случаях с оксидом и без него (рис. 7).

Рис. 7. Распределение полей отраженного излучения при А/Л to 1.25 (верхняя строка) и А/Л to 1.0 (нижняя строка) для решетки без слоя оксида (левый столбец) и с ним (правый столбец)

Согласно зависимостям на рис. 6 основной вклад в «равномерное» распределение поля (как в левом столбце на рис. 7) дает нулевая мода, а неоднородности обеспечиваются модами более высоких порядков, в частности ±1-ми, (правый график на рис. 6).

При фиксированном периоде (Л ~ Л/1.25 = 1.24 мкм) построены зависимости коэффициентов отражения от высоты решетки (при L_ox = 0.2 мкм) и толщины оксида (при Һ = 50 нм) (рис. 8).

Рис. 8. Зависимости модулей коэффициентов отражения |^_м0| от высотьi решетки Һ при L_Ox = 0.2 мкм (слева) и от толщины оксида L_ox при Һ = 50 нм (справа). Штриховыми линиями указаны точки для построения распределений поля

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

Рис. 9. Распределение полей отраженного излучения при высоте решетки Һ = 50, 150 нм и толщине оксида L_ox = 0.2, 1.0 мкм. Обратим еще раз внимание на вклад однородной нулевой моды и неоднородных мод больших порядков в «равномерность» распределения

• 4.    Заключение

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